Line data Source code
1 : /* Copyright (C) 2000 The PARI group.
2 :
3 : This file is part of the PARI/GP package.
4 :
5 : PARI/GP is free software; you can redistribute it and/or modify it under the
6 : terms of the GNU General Public License as published by the Free Software
7 : Foundation. It is distributed in the hope that it will be useful, but WITHOUT
8 : ANY WARRANTY WHATSOEVER.
9 :
10 : Check the License for details. You should have received a copy of it, along
11 : with the package; see the file 'COPYING'. If not, write to the Free Software
12 : Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. */
13 :
14 : /***********************************************************************/
15 : /** **/
16 : /** ARITHMETIC OPERATIONS ON POLYNOMIALS **/
17 : /** (second part) **/
18 : /** **/
19 : /***********************************************************************/
20 : #include "pari.h"
21 : #include "paripriv.h"
22 :
23 : /* compute Newton sums S_1(P), ... , S_n(P). S_k(P) = sum a_j^k, a_j root of P
24 : * If N != NULL, assume p-adic roots and compute mod N [assume integer coeffs]
25 : * If T != NULL, compute mod (T,N) [assume integer coeffs if N != NULL]
26 : * If y0!= NULL, precomputed i-th powers, i=1..m, m = length(y0).
27 : * Not memory clean in the latter case */
28 : GEN
29 29540 : polsym_gen(GEN P, GEN y0, long n, GEN T, GEN N)
30 : {
31 29540 : long dP=degpol(P), i, k, m;
32 : pari_sp av1, av2;
33 : GEN s,y,P_lead;
34 :
35 29540 : if (n<0) pari_err_IMPL("polsym of a negative n");
36 29540 : if (typ(P) != t_POL) pari_err_TYPE("polsym",P);
37 29540 : if (!signe(P)) pari_err_ROOTS0("polsym");
38 29540 : y = cgetg(n+2,t_COL);
39 29540 : if (y0)
40 : {
41 9107 : if (typ(y0) != t_COL) pari_err_TYPE("polsym_gen",y0);
42 9107 : m = lg(y0)-1;
43 9107 : for (i=1; i<=m; i++) gel(y,i) = gel(y0,i); /* not memory clean */
44 : }
45 : else
46 : {
47 20433 : m = 1;
48 20433 : gel(y,1) = stoi(dP);
49 : }
50 29540 : P += 2; /* strip codewords */
51 :
52 29540 : P_lead = gel(P,dP); if (gequal1(P_lead)) P_lead = NULL;
53 29540 : if (P_lead)
54 : {
55 7 : if (N) P_lead = Fq_inv(P_lead,T,N);
56 7 : else if (T) P_lead = QXQ_inv(P_lead,T);
57 : }
58 78211 : for (k=m; k<=n; k++)
59 : {
60 48671 : av1 = avma; s = (dP>=k)? gmulsg(k,gel(P,dP-k)): gen_0;
61 216762 : for (i=1; i<k && i<=dP; i++)
62 168091 : s = gadd(s, gmul(gel(y,k-i+1),gel(P,dP-i)));
63 48671 : if (N)
64 : {
65 13412 : s = Fq_red(s, T, N);
66 13412 : if (P_lead) s = Fq_mul(s, P_lead, T, N);
67 : }
68 35259 : else if (T)
69 : {
70 0 : s = grem(s, T);
71 0 : if (P_lead) s = grem(gmul(s, P_lead), T);
72 : }
73 : else
74 35259 : if (P_lead) s = gdiv(s, P_lead);
75 48671 : av2 = avma; gel(y,k+1) = gerepile(av1,av2, gneg(s));
76 : }
77 29540 : return y;
78 : }
79 :
80 : GEN
81 17472 : polsym(GEN x, long n)
82 : {
83 17472 : return polsym_gen(x, NULL, n, NULL,NULL);
84 : }
85 :
86 : /* centered residue x mod p. po2 = shifti(p, -1) or NULL (euclidean residue) */
87 : GEN
88 86341568 : centermodii(GEN x, GEN p, GEN po2)
89 : {
90 86341568 : GEN y = remii(x, p);
91 86346745 : switch(signe(y))
92 : {
93 13778714 : case 0: break;
94 47618011 : case 1: if (po2 && abscmpii(y,po2) > 0) y = subii(y, p);
95 47553748 : break;
96 25092523 : case -1: if (!po2 || abscmpii(y,po2) > 0) y = addii(y, p);
97 25057699 : break;
98 : }
99 86247658 : return y;
100 : }
101 :
102 : static long
103 0 : s_centermod(long x, ulong pp, ulong pps2)
104 : {
105 0 : long y = x % (long)pp;
106 0 : if (y < 0) y += pp;
107 0 : return Fl_center(y, pp,pps2);
108 : }
109 :
110 : /* for internal use */
111 : GEN
112 8368338 : centermod_i(GEN x, GEN p, GEN ps2)
113 : {
114 : long i, lx;
115 : pari_sp av;
116 : GEN y;
117 :
118 8368338 : if (!ps2) ps2 = shifti(p,-1);
119 8368717 : switch(typ(x))
120 : {
121 3670688 : case t_INT: return centermodii(x,p,ps2);
122 :
123 3897612 : case t_POL: lx = lg(x);
124 3897612 : y = cgetg(lx,t_POL); y[1] = x[1];
125 32062056 : for (i=2; i<lx; i++)
126 : {
127 28163639 : av = avma;
128 28163639 : gel(y,i) = gerepileuptoint(av, centermodii(gel(x,i),p,ps2));
129 : }
130 3898417 : return normalizepol_lg(y, lx);
131 :
132 799577 : case t_COL: lx = lg(x);
133 799577 : y = cgetg(lx,t_COL);
134 799578 : for (i=1; i<lx; i++) gel(y,i) = centermodii(gel(x,i),p,ps2);
135 799575 : return y;
136 :
137 840 : case t_MAT: lx = lg(x);
138 840 : y = cgetg(lx,t_MAT);
139 840 : for (i=1; i<lx; i++) gel(y,i) = centermod_i(gel(x,i),p,ps2);
140 840 : return y;
141 :
142 0 : case t_VECSMALL: lx = lg(x);
143 : {
144 0 : ulong pp = itou(p), pps2 = itou(ps2);
145 0 : y = cgetg(lx,t_VECSMALL);
146 0 : for (i=1; i<lx; i++) y[i] = s_centermod(x[i], pp, pps2);
147 0 : return y;
148 : }
149 : }
150 0 : return x;
151 : }
152 :
153 : GEN
154 4802052 : centermod(GEN x, GEN p) { return centermod_i(x,p,NULL); }
155 :
156 : /***********************************************************************/
157 : /** **/
158 : /** FACTORIZATION **/
159 : /** **/
160 : /***********************************************************************/
161 : #define assign_or_fail(x,y) { GEN __x = x;\
162 : if (!*y) *y=__x; else if (!gequal(__x,*y)) return 0;\
163 : }
164 : #define update_prec(x,y) { long __x = x; if (__x < *y) *y=__x; }
165 :
166 : static const long tsh = 6;
167 : static long
168 14588 : code(long t1, long t2) { return (t1 << tsh) | t2; }
169 : void
170 14826 : RgX_type_decode(long x, long *t1, long *t2)
171 : {
172 14826 : *t1 = x >> tsh;
173 14826 : *t2 = (x & ((1L<<tsh)-1));
174 14826 : }
175 : int
176 804385 : RgX_type_is_composite(long t) { return t >= tsh; }
177 :
178 : static int
179 3107848 : settype(GEN c, long *t, GEN *p, GEN *pol, long *pa, GEN *ff, long *t2)
180 : {
181 : long j;
182 3107848 : switch(typ(c))
183 : {
184 : case t_INT: case t_FRAC:
185 2075353 : break;
186 : case t_REAL:
187 7499 : update_prec(precision(c), pa);
188 7499 : t[2]=1; break;
189 : case t_INTMOD:
190 1519 : assign_or_fail(gel(c,1),p);
191 1519 : t[3]=1; break;
192 : case t_FFELT:
193 462959 : if (!*ff) *ff=c; else if (!FF_samefield(c,*ff)) return 0;
194 462959 : assign_or_fail(FF_p_i(c),p);
195 462959 : t[5]=1; break;
196 : case t_COMPLEX:
197 2800 : for (j=1; j<=2; j++)
198 : {
199 1869 : GEN d = gel(c,j);
200 1869 : switch(typ(d))
201 : {
202 : case t_INT: case t_FRAC:
203 1302 : t[4]=1; break;
204 : case t_REAL:
205 546 : update_prec(precision(d), pa);
206 546 : t[6]=1; break;
207 : case t_INTMOD:
208 14 : assign_or_fail(gel(d,1),p);
209 14 : if (!signe(*p) || mod4(*p) != 3) return 0;
210 7 : if (!*t2) *t2 = t_COMPLEX;
211 7 : t[3]=1; break;
212 : case t_PADIC:
213 7 : update_prec(precp(d)+valp(d), pa);
214 7 : assign_or_fail(gel(d,2),p);
215 7 : if (!*t2) *t2 = t_COMPLEX;
216 7 : t[7]=1; break;
217 0 : default: return 0;
218 : }
219 : }
220 931 : if (!t[6]) assign_or_fail(mkpoln(3, gen_1,gen_0,gen_1), pol); /*x^2+1*/
221 931 : break;
222 : case t_PADIC:
223 77 : update_prec(precp(c)+valp(c), pa);
224 77 : assign_or_fail(gel(c,2),p);
225 77 : t[7]=1; break;
226 : case t_QUAD:
227 28 : assign_or_fail(gel(c,1),pol);
228 84 : for (j=2; j<=3; j++)
229 : {
230 56 : GEN d = gel(c,j);
231 56 : switch(typ(d))
232 : {
233 : case t_INT: case t_FRAC:
234 21 : t[8]=1; break;
235 : case t_INTMOD:
236 28 : assign_or_fail(gel(d,1),p);
237 28 : if (*t2 != t_POLMOD) *t2 = t_QUAD;
238 28 : t[3]=1; break;
239 : case t_PADIC:
240 7 : update_prec(precp(d)+valp(d), pa);
241 7 : assign_or_fail(gel(d,2),p);
242 7 : if (*t2 != t_POLMOD) *t2 = t_QUAD;
243 7 : t[7]=1; break;
244 0 : default: return 0;
245 : }
246 : }
247 28 : break;
248 : case t_POLMOD:
249 75124 : assign_or_fail(gel(c,1),pol);
250 450744 : for (j=1; j<=2; j++)
251 : {
252 150248 : GEN pbis = NULL, polbis = NULL;
253 : long pabis;
254 150248 : switch(RgX_type(gel(c,j),&pbis,&polbis,&pabis))
255 : {
256 150227 : case t_INT: t[9]=1; break;
257 14 : case t_INTMOD: t[3]=1; *t2 = t_POLMOD; break;
258 7 : case t_PADIC: t[7]=1; *t2 = t_POLMOD; update_prec(pabis,pa); break;
259 0 : default: return 0;
260 : }
261 150248 : if (pbis) assign_or_fail(pbis,p);
262 150248 : if (polbis) assign_or_fail(polbis,pol);
263 : }
264 75124 : break;
265 484351 : default: return 0;
266 : }
267 2623490 : return 1;
268 : }
269 : /* t[0..1] unused. Other values, if set, indicate a coefficient of type
270 : * t[2] : t_REAL
271 : * t[3] : t_INTMOD
272 : * t[4] : t_COMPLEX of rationals (t_INT/t_FRAC)
273 : * t[5] : t_FFELT
274 : * t[6] : t_COMPLEX of t_REAL
275 : * t[7] : t_PADIC
276 : * t[8] : t_QUAD of rationals (t_INT/t_FRAC)
277 : * t[9]: t_POLMOD of rationals (t_INT/t_FRAC) */
278 : /* if t2 != 0: t_POLMOD/t_QUAD/t_COMPLEX of modular (t_INTMOD/t_PADIC,
279 : * given by t) */
280 : static long
281 575527 : choosetype(long *t, long t2, GEN ff, GEN *pol)
282 : {
283 575527 : if (t[5]) /* ffelt */
284 : {
285 103705 : if (t2 ||t[2]||t[4]||t[6]||t[8]||t[9]) return 0;
286 103705 : *pol=ff; return t_FFELT;
287 : }
288 471822 : if (t[6]) /* inexact, complex */
289 : {
290 280 : if (t2 ||t[3]||t[7]||t[9]) return 0;
291 280 : return t_COMPLEX;
292 : }
293 471542 : if (t[2]) /* inexact, real */
294 : {
295 5833 : if (t2 ||t[3]||t[7]||t[9]) return 0;
296 5833 : return t[4]?t_COMPLEX:t_REAL;
297 : }
298 465709 : if (t2) /* polmod/quad/complex of intmod/padic */
299 : {
300 56 : if (t[3]) return code(t2,t_INTMOD);
301 21 : if (t[7]) return code(t2,t_PADIC);
302 : }
303 465653 : if (t[9]) return code(t_POLMOD,t_INT);
304 451730 : if (t[8]) return code(t_QUAD,t_INT);
305 451723 : if (t[4]) return code(t_COMPLEX,t_INT);
306 451121 : if (t[3]) return t_INTMOD;
307 450750 : if (t[7]) return t_PADIC;
308 450715 : return t_INT;
309 : }
310 : long
311 1059948 : RgX_type(GEN x, GEN *p, GEN *pol, long *pa)
312 : {
313 1059948 : long t[] = {0,0,0,0,0,0,0,0,0,0};
314 1059948 : long tx = typ(x), t2 = 0;
315 1059948 : GEN ff = NULL;
316 1059948 : *p = *pol = NULL; *pa = LONG_MAX;
317 1059948 : if (is_scalar_t(tx))
318 : {
319 10096 : if (tx == t_POLMOD) return 0;
320 10033 : if (!settype(x,t,p,pol,pa,&ff,&t2)) return 0;
321 : }
322 1049852 : else if (tx == t_MAT)
323 : {
324 7 : long j, lx = lg(x);
325 21 : for (j = 1; j < lx; j++)
326 : {
327 14 : GEN c = gel(x,j);
328 14 : long i, l = lg(c);
329 42 : for (i=1; i<l; i++)
330 28 : if (!settype(gel(c,i),t,p,pol,pa,&ff,&t2)) return 0;
331 : }
332 : }
333 : else /* t_POL, t_SER */
334 : {
335 1049845 : long i, lx = lg(x);
336 3663274 : for (i=2; i<lx; i++)
337 3097787 : if (!settype(gel(x,i),t,p,pol,pa,&ff,&t2)) return 0;
338 : }
339 575527 : return choosetype(t,t2,ff,pol);
340 : }
341 :
342 : GEN
343 0 : factor0(GEN x,long flag)
344 : {
345 0 : return (flag<0)? factor(x): boundfact(x,flag);
346 : }
347 :
348 : /* only present for interface with GP */
349 : GEN
350 37017 : gp_factor0(GEN x, GEN flag)
351 : {
352 : ulong B;
353 37017 : if (!flag) return factor(x);
354 35 : if (typ(flag) != t_INT || signe(flag) < 0) pari_err_FLAG("factor");
355 35 : switch(lgefint(flag))
356 : {
357 7 : case 2: B = 0; break;
358 28 : case 3: B = flag[2]; break;
359 0 : default: pari_err_OVERFLOW("factor [large prime bound]");
360 : return NULL; /*LCOV_EXCL_LINE*/
361 : }
362 35 : return boundfact(x, B);
363 : }
364 :
365 : GEN
366 82709 : deg1_from_roots(GEN L, long v)
367 : {
368 82709 : long i, l = lg(L);
369 82709 : GEN z = cgetg(l,t_COL);
370 188894 : for (i=1; i<l; i++)
371 106185 : gel(z,i) = deg1pol_shallow(gen_1, gneg(gel(L,i)), v);
372 82709 : return z;
373 : }
374 : GEN
375 980 : roots_from_deg1(GEN x)
376 : {
377 980 : long i,l = lg(x);
378 980 : GEN r = cgetg(l,t_VEC);
379 980 : for (i=1; i<l; i++) { GEN P = gel(x,i); gel(r,i) = gneg(gel(P,2)); }
380 980 : return r;
381 : }
382 :
383 : static GEN
384 35 : gauss_factor_p(GEN p)
385 : {
386 35 : GEN a, b; (void)cornacchia(gen_1, p, &a,&b);
387 35 : return mkcomplex(a, b);
388 : }
389 :
390 : static GEN
391 42 : gauss_primpart(GEN x, GEN *c)
392 : {
393 42 : GEN a = gel(x,1), b = gel(x,2), n = gcdii(a, b);
394 42 : *c = n; if (n == gen_1) return x;
395 42 : retmkcomplex(diviiexact(a,n), diviiexact(b,n));
396 : }
397 :
398 : static GEN
399 70 : gauss_primpart_try(GEN x, GEN c)
400 : {
401 : GEN r, y;
402 70 : if (typ(x) == t_INT)
403 : {
404 42 : y = dvmdii(x, c, &r); if (r != gen_0) return NULL;
405 : }
406 : else
407 : {
408 28 : GEN a = gel(x,1), b = gel(x,2); y = cgetg(3, t_COMPLEX);
409 28 : gel(y,1) = dvmdii(a, c, &r); if (r != gen_0) return NULL;
410 14 : gel(y,2) = dvmdii(b, c, &r); if (r != gen_0) return NULL;
411 : }
412 56 : return y;
413 : }
414 :
415 : static int
416 77 : gauss_cmp(GEN x, GEN y)
417 : {
418 : int v;
419 77 : if (typ(x) != t_COMPLEX)
420 0 : return (typ(y) == t_COMPLEX)? -1: gcmp(x, y);
421 77 : if (typ(y) != t_COMPLEX) return 1;
422 49 : v = cmpii(gel(x,2), gel(y,2));
423 49 : if (v) return v;
424 21 : return gcmp(gel(x,1), gel(y,1));
425 : }
426 :
427 : /* 0 or canonical representative in Z[i]^* / <i> (impose imag(x) >= 0) */
428 : static GEN
429 91 : gauss_normal(GEN x)
430 : {
431 91 : if (typ(x) != t_COMPLEX) return (signe(x) < 0)? absi(x): x;
432 84 : if (signe(gel(x,1)) < 0) x = gneg(x);
433 84 : if (signe(gel(x,2)) < 0) x = mulcxI(x);
434 84 : return x;
435 : }
436 :
437 : static GEN
438 42 : gauss_factor(GEN x)
439 : {
440 42 : pari_sp av = avma;
441 42 : GEN a = gel(x,1), b = gel(x,2), d = gen_1, n, y, fa, P, E, P2, E2;
442 42 : long t1 = typ(a);
443 42 : long t2 = typ(b), i, j, l, exp = 0;
444 42 : if (t1 == t_FRAC) d = gel(a,2);
445 42 : if (t2 == t_FRAC) d = lcmii(d, gel(b,2));
446 42 : if (d == gen_1) y = x;
447 : else
448 : {
449 14 : y = gmul(x, d);
450 14 : a = gel(y,1); t1 = typ(a);
451 14 : b = gel(y,2); t2 = typ(b);
452 : }
453 42 : if (t1 != t_INT || t2 != t_INT) return NULL;
454 42 : y = gauss_primpart(y, &n);
455 42 : fa = factor(cxnorm(y));
456 42 : P = gel(fa,1);
457 42 : E = gel(fa,2); l = lg(P);
458 42 : P2 = cgetg(l, t_COL);
459 42 : E2 = cgetg(l, t_COL);
460 98 : for (j = 1, i = l-1; i > 0; i--) /* remove largest factors first */
461 : { /* either p = 2 (ramified) or those factors split in Q(i) */
462 56 : GEN p = gel(P,i), w, w2, t, we, pe;
463 56 : long v, e = itos(gel(E,i));
464 56 : int is2 = absequaliu(p, 2);
465 56 : w = is2? mkcomplex(gen_1,gen_1): gauss_factor_p(p);
466 56 : w2 = gauss_normal( gconj(w) );
467 : /* w * w2 * I^3 = p, w2 = gconj(w) * I */
468 56 : pe = powiu(p, e);
469 56 : we = gpowgs(w, e);
470 56 : t = gauss_primpart_try( gmul(y, gconj(we)), pe );
471 56 : if (t) y = t; /* y /= w^e */
472 : else {
473 : /* y /= conj(w)^e, should be y /= w2^e */
474 14 : y = gauss_primpart_try( gmul(y, we), pe );
475 14 : swap(w, w2); exp -= e; /* += 3*e mod 4 */
476 : }
477 56 : gel(P,i) = w;
478 56 : v = Z_pvalrem(n, p, &n);
479 56 : if (v) {
480 7 : exp -= v; /* += 3*v mod 4 */
481 7 : if (is2) v <<= 1; /* 2 = w^2 I^3 */
482 : else {
483 0 : gel(P2,j) = w2;
484 0 : gel(E2,j) = utoipos(v); j++;
485 : }
486 7 : gel(E,i) = stoi(e + v);
487 : }
488 56 : v = Z_pvalrem(d, p, &d);
489 56 : if (v) {
490 7 : exp += v; /* -= 3*v mod 4 */
491 7 : if (is2) v <<= 1; /* 2 is ramified */
492 : else {
493 7 : gel(P2,j) = w2;
494 7 : gel(E2,j) = utoineg(v); j++;
495 : }
496 7 : gel(E,i) = stoi(e - v);
497 : }
498 56 : exp &= 3;
499 : }
500 42 : if (j > 1) {
501 7 : long k = 1;
502 7 : GEN P1 = cgetg(l, t_COL);
503 7 : GEN E1 = cgetg(l, t_COL);
504 : /* remove factors with exponent 0 */
505 14 : for (i = 1; i < l; i++)
506 7 : if (signe(gel(E,i)))
507 : {
508 0 : gel(P1,k) = gel(P,i);
509 0 : gel(E1,k) = gel(E,i);
510 0 : k++;
511 : }
512 7 : setlg(P1, k); setlg(E1, k);
513 7 : setlg(P2, j); setlg(E2, j);
514 7 : fa = famat_mul_shallow(mkmat2(P1,E1), mkmat2(P2,E2));
515 : }
516 42 : if (!is_pm1(n) || !is_pm1(d))
517 : {
518 21 : GEN Fa = factor(gdiv(n, d));
519 21 : P = gel(Fa,1); l = lg(P);
520 21 : E = gel(Fa,2);
521 49 : for (i = 1; i < l; i++)
522 : {
523 28 : GEN w, p = gel(P,i);
524 : long e;
525 : int is2;
526 28 : switch(mod4(p))
527 : {
528 14 : case 3: continue;
529 7 : case 2: is2 = 1; break;
530 7 : default:is2 = 0; break;
531 : }
532 14 : e = itos(gel(E,i));
533 14 : w = is2? mkcomplex(gen_1,gen_1): gauss_factor_p(p);
534 14 : gel(P,i) = w;
535 14 : if (is2)
536 7 : gel(E,i) = stoi(2*e);
537 : else
538 : {
539 7 : P = shallowconcat(P, gauss_normal( gconj(w) ));
540 7 : E = shallowconcat(E, gel(E,i));
541 : }
542 14 : exp -= e; /* += 3*e mod 4 */
543 14 : exp &= 3;
544 : }
545 21 : gel(Fa,1) = P;
546 21 : gel(Fa,2) = E;
547 21 : fa = famat_mul_shallow(fa, Fa);
548 : }
549 42 : fa = sort_factor(fa, (void*)&gauss_cmp, &cmp_nodata);
550 :
551 42 : y = gmul(y, powIs(exp));
552 42 : if (!gequal1(y)) {
553 35 : gel(fa,1) = shallowconcat(mkcol(y), gel(fa,1));
554 35 : gel(fa,2) = shallowconcat(gen_1, gel(fa,2));
555 : }
556 42 : return gerepilecopy(av, fa);
557 : }
558 :
559 : GEN
560 41182 : factor(GEN x)
561 : {
562 41182 : long tx=typ(x), lx, i, pa, v, r1;
563 : pari_sp av, tetpil;
564 : GEN y, p, p1, p2, pol;
565 :
566 41182 : if (gequal0(x))
567 56 : switch(tx)
568 : {
569 : case t_INT:
570 : case t_COMPLEX:
571 : case t_POL:
572 56 : case t_RFRAC: return prime_fact(x);
573 0 : default: pari_err_TYPE("factor",x);
574 : }
575 41126 : av = avma;
576 41126 : switch(tx)
577 : {
578 39047 : case t_INT: return Z_factor(x);
579 :
580 : case t_FRAC:
581 182 : p1 = Z_factor(gel(x,1));
582 182 : p2 = Z_factor(gel(x,2)); gel(p2,2) = ZC_neg(gel(p2,2));
583 182 : return gerepilecopy(av, merge_factor(p1,p2,(void*)&cmpii,cmp_nodata));
584 :
585 : case t_POL:
586 1848 : tx=RgX_type(x,&p,&pol,&pa);
587 1848 : switch(tx)
588 : {
589 7 : case 0: pari_err_IMPL("factor for general polynomials");
590 1610 : case t_INT: return QX_factor(x);
591 7 : case t_INTMOD: return factmod(x,p);
592 :
593 7 : case t_COMPLEX: y=cgetg(3,t_MAT); lx=lg(x)-2;
594 7 : av = avma; p1 = roots(x,pa); tetpil = avma;
595 7 : p1 = deg1_from_roots(p1, varn(x));
596 7 : gel(y,1) = gerepile(av,tetpil,p1);
597 7 : gel(y,2) = const_col(lx-1, gen_1); return y;
598 :
599 14 : case t_REAL: y=cgetg(3,t_MAT); lx=lg(x)-2; v=varn(x);
600 14 : av=avma; p1=cleanroots(x,pa); tetpil=avma;
601 35 : for(r1=1; r1<lx; r1++)
602 28 : if (typ(gel(p1,r1)) == t_COMPLEX) break;
603 14 : lx=(r1+lx)>>1; p2=cgetg(lx,t_COL);
604 35 : for(i=1; i<r1; i++)
605 21 : gel(p2,i) = deg1pol_shallow(gen_1, negr(gel(p1,i)), v);
606 21 : for( ; i<lx; i++)
607 : {
608 7 : GEN a = gel(p1,2*i-r1);
609 7 : p = cgetg(5, t_POL); gel(p2,i) = p;
610 7 : p[1] = x[1];
611 7 : gel(p,2) = gnorm(a);
612 7 : gel(p,3) = gmul2n(gel(a,1),1); togglesign(gel(p,3));
613 7 : gel(p,4) = gen_1;
614 : }
615 14 : gel(y,1) = gerepile(av,tetpil,p2);
616 14 : gel(y,2) = const_col(lx-1, gen_1); return y;
617 :
618 21 : case t_PADIC: return factorpadic(x,p,pa);
619 :
620 56 : case t_FFELT: return FFX_factor(x,pol);
621 :
622 : default:
623 : {
624 : GEN w;
625 126 : long killv, t1, t2, v2 = varn(pol);
626 126 : x = leafcopy(x); lx=lg(x);
627 126 : pol = leafcopy(pol);
628 126 : v = pari_var_next_temp();
629 1001 : for(i=2; i<lx; i++)
630 : {
631 875 : p1 = gel(x,i);
632 875 : switch(typ(p1))
633 : {
634 28 : case t_QUAD: p1++;
635 : case t_COMPLEX:
636 49 : gel(x,i) = mkpolmod(deg1pol_shallow(gel(p1,2), gel(p1,1), v), pol);
637 : }
638 : }
639 126 : killv = (avma != (pari_sp)pol);
640 126 : if (killv) setvarn(pol, fetch_var());
641 126 : RgX_type_decode(tx, &t1, &t2);
642 126 : switch (t2)
643 : {
644 70 : case t_INT: p1 = nffactor(pol,x); break;
645 : case t_INTMOD:
646 35 : pol = RgX_to_FpX(pol, p);
647 35 : if (FpX_is_squarefree(pol,p) && FpX_nbfact(pol, p) == 1)
648 : {
649 21 : p1 = factorff(x,p,pol); break;
650 : }
651 : /*fall through*/
652 35 : default: pari_err_IMPL("factor for general polynomial");
653 : return NULL; /* LCOV_EXCL_LINE */
654 : }
655 91 : if (killv) setvarn(pol, v2);
656 91 : switch (t1)
657 : {
658 : case t_POLMOD:
659 63 : if (killv) (void)delete_var();
660 63 : return gerepileupto(av,p1);
661 14 : case t_COMPLEX: w = gen_I(); break;
662 14 : case t_QUAD: w = mkquad(pol,gen_0,gen_1);
663 14 : break;
664 0 : default: pari_err_IMPL("factor for general polynomial");
665 : return NULL; /* LCOV_EXCL_LINE */
666 : }
667 28 : p2=gel(p1,1);
668 70 : for(i=1; i<lg(p2); i++)
669 : {
670 42 : GEN P = gel(p2,i);
671 : long j;
672 140 : for(j=2; j<lg(P); j++)
673 : {
674 98 : GEN p = gel(P,j);
675 98 : if(typ(p)==t_POLMOD) gel(P,j) = gsubst(gel(p,2),v,w);
676 : }
677 : }
678 28 : if (killv) (void)delete_var();
679 28 : return gerepilecopy(av, p1);
680 : }
681 : }
682 : case t_RFRAC: {
683 7 : GEN a = gel(x,1), b = gel(x,2);
684 7 : y = factor(b); gel(y,2) = gneg_i(gel(y,2));
685 7 : if (typ(a)==t_POL && varn(a)==varn(b)) y = famat_mul_shallow(factor(a), y);
686 7 : return gerepilecopy(av, y);
687 : }
688 :
689 : case t_COMPLEX:
690 42 : y = gauss_factor(x);
691 42 : if (y) return y;
692 : /* fall through */
693 : }
694 0 : pari_err_TYPE("factor",x);
695 : return NULL; /* LCOV_EXCL_LINE */
696 : }
697 :
698 : /*******************************************************************/
699 : /* */
700 : /* ROOTS --> MONIC POLYNOMIAL */
701 : /* */
702 : /*******************************************************************/
703 : static GEN
704 414604 : normalized_mul(void *E, GEN x, GEN y)
705 : {
706 414604 : long a = gel(x,1)[1], b = gel(y,1)[1];
707 : (void) E;
708 829208 : return mkvec2(mkvecsmall(a + b),
709 829208 : RgX_mul_normalized(gel(x,2),a, gel(y,2),b));
710 : }
711 : /* L = [Vecsmall([a]), A], with a > 0, A an RgX, deg(A) < a; return X^a + A */
712 : static GEN
713 264723 : normalized_to_RgX(GEN L)
714 : {
715 264723 : long i, a = gel(L,1)[1];
716 264723 : GEN A = gel(L,2);
717 264723 : GEN z = cgetg(a + 3, t_POL);
718 264723 : z[1] = evalsigne(1) | evalvarn(varn(A));
719 264723 : for (i = 2; i < lg(A); i++) gel(z,i) = gcopy(gel(A,i));
720 264723 : for ( ; i < a+2; i++) gel(z,i) = gen_0;
721 264723 : gel(z,i) = gen_1; return z;
722 : }
723 :
724 : /* compute prod (x - a[i]) */
725 : GEN
726 239511 : roots_to_pol(GEN a, long v)
727 : {
728 239511 : pari_sp av = avma;
729 239511 : long i, k, lx = lg(a);
730 : GEN L;
731 239511 : if (lx == 1) return pol_1(v);
732 239476 : L = cgetg(lx, t_VEC);
733 512043 : for (k=1,i=1; i<lx-1; i+=2)
734 : {
735 272567 : GEN s = gel(a,i), t = gel(a,i+1);
736 272567 : GEN x0 = gmul(s,t);
737 272567 : GEN x1 = gneg(gadd(s,t));
738 272567 : gel(L,k++) = mkvec2(mkvecsmall(2), deg1pol_shallow(x1,x0,v));
739 : }
740 466202 : if (i < lx) gel(L,k++) = mkvec2(mkvecsmall(1),
741 226726 : scalarpol_shallow(gneg(gel(a,i)), v));
742 239476 : setlg(L, k); L = gen_product(L, NULL, normalized_mul);
743 239476 : return gerepileupto(av, normalized_to_RgX(L));
744 : }
745 :
746 : /* prod_{i=1..r1} (x - a[i]) prod_{i=1..r2} (x - a[i])(x - conj(a[i]))*/
747 : GEN
748 25247 : roots_to_pol_r1(GEN a, long v, long r1)
749 : {
750 25247 : pari_sp av = avma;
751 25247 : long i, k, lx = lg(a);
752 : GEN L;
753 25247 : if (lx == 1) return pol_1(v);
754 25247 : L = cgetg(lx, t_VEC);
755 123158 : for (k=1,i=1; i<r1; i+=2)
756 : {
757 97911 : GEN s = gel(a,i), t = gel(a,i+1);
758 97911 : GEN x0 = gmul(s,t);
759 97911 : GEN x1 = gneg(gadd(s,t));
760 97911 : gel(L,k++) = mkvec2(mkvecsmall(2), deg1pol_shallow(x1,x0,v));
761 : }
762 31253 : if (i < r1+1) gel(L,k++) = mkvec2(mkvecsmall(1),
763 6006 : scalarpol_shallow(gneg(gel(a,i)), v));
764 101364 : for (i=r1+1; i<lx; i++)
765 : {
766 76117 : GEN s = gel(a,i);
767 76117 : GEN x0 = gnorm(s);
768 76117 : GEN x1 = gneg(gtrace(s));
769 76117 : gel(L,k++) = mkvec2(mkvecsmall(2), deg1pol_shallow(x1,x0,v));
770 : }
771 25247 : setlg(L, k); L = gen_product(L, NULL, normalized_mul);
772 25247 : return gerepileupto(av, normalized_to_RgX(L));
773 : }
774 :
775 : /*******************************************************************/
776 : /* */
777 : /* FACTORBACK */
778 : /* */
779 : /*******************************************************************/
780 : static GEN
781 21 : idmulred(void *nf, GEN x, GEN y) { return idealmulred((GEN) nf, x, y); }
782 : static GEN
783 623 : idpowred(void *nf, GEN x, GEN n) { return idealpowred((GEN) nf, x, n); }
784 : static GEN
785 2142 : idmul(void *nf, GEN x, GEN y) { return idealmul((GEN) nf, x, y); }
786 : static GEN
787 3283 : idpow(void *nf, GEN x, GEN n) { return idealpow((GEN) nf, x, n); }
788 : static GEN
789 4165 : eltmul(void *nf, GEN x, GEN y) { return nfmul((GEN) nf, x, y); }
790 : static GEN
791 16962 : eltpow(void *nf, GEN x, GEN n) { return nfpow((GEN) nf, x, n); }
792 : static GEN
793 35637233 : mul(void *a, GEN x, GEN y) { (void)a; return gmul(x,y);}
794 : static GEN
795 51280993 : powi(void *a, GEN x, GEN y) { (void)a; return powgi(x,y);}
796 : static GEN
797 23870 : Fpmul(void *a, GEN x, GEN y) { return Fp_mul(x,y,(GEN)a); }
798 : static GEN
799 120190 : Fppow(void *a, GEN x, GEN n) { return Fp_pow(x,n,(GEN)a); }
800 :
801 : /* [L,e] = [fa, NULL] or [elts, NULL] or [elts, exponents] */
802 : GEN
803 15791817 : gen_factorback(GEN L, GEN e, GEN (*_mul)(void*,GEN,GEN),
804 : GEN (*_pow)(void*,GEN,GEN), void *data)
805 : {
806 15791817 : pari_sp av = avma;
807 : long k, l, lx;
808 : GEN p,x;
809 :
810 15791817 : if (e) /* supplied vector of exponents */
811 600761 : p = L;
812 : else
813 : {
814 15191056 : switch(typ(L)) {
815 : case t_VEC:
816 : case t_COL: /* product of the L[i] */
817 244 : return gerepileupto(av, gen_product(L, data, _mul));
818 : case t_MAT: /* genuine factorization */
819 15190812 : l = lg(L);
820 15190812 : if (l == 1) return gen_1;
821 15190805 : if (l == 3) break;
822 : /*fall through*/
823 : default:
824 7 : pari_err_TYPE("factorback [not a factorization]", L);
825 : }
826 15190798 : p = gel(L,1);
827 15190798 : e = gel(L,2);
828 : }
829 : /* p = elts, e = expo */
830 15791559 : lx = lg(p);
831 : /* check whether e is an integral vector of correct length */
832 15791559 : switch(typ(e))
833 : {
834 : case t_VECSMALL:
835 252 : if (lx != lg(e))
836 0 : pari_err_TYPE("factorback [not an exponent vector]", e);
837 252 : if (lx == 1) return gen_1;
838 252 : x = cgetg(lx,t_VEC);
839 854 : for (l=1,k=1; k<lx; k++)
840 602 : if (e[k]) gel(x,l++) = _pow(data, gel(p,k), stoi(e[k]));
841 252 : break;
842 : case t_VEC: case t_COL:
843 15791307 : if (lx != lg(e) || !RgV_is_ZV(e))
844 7 : pari_err_TYPE("factorback [not an exponent vector]", e);
845 15791300 : if (lx == 1) return gen_1;
846 15780000 : x = cgetg(lx,t_VEC);
847 67364591 : for (l=1,k=1; k<lx; k++)
848 51584591 : if (signe(gel(e,k))) gel(x,l++) = _pow(data, gel(p,k), gel(e,k));
849 15780000 : break;
850 : default:
851 0 : pari_err_TYPE("factorback [not an exponent vector]", e);
852 0 : return NULL;
853 : }
854 15780252 : x[0] = evaltyp(t_VEC) | _evallg(l);
855 15780252 : return gerepileupto(av, gen_product(x, data, _mul));
856 : }
857 :
858 : GEN
859 1799 : idealfactorback(GEN nf, GEN L, GEN e, int red)
860 : {
861 1799 : nf = checknf(nf);
862 1799 : if (red) return gen_factorback(L, e, &idmulred, &idpowred, (void*)nf);
863 1197 : else return gen_factorback(L, e, &idmul, &idpow, (void*)nf);
864 : }
865 :
866 : GEN
867 13307 : nffactorback(GEN nf, GEN L, GEN e)
868 13307 : { return gen_factorback(L, e, &eltmul, &eltpow, (void*)checknf(nf)); }
869 :
870 : GEN
871 96320 : FpV_factorback(GEN L, GEN e, GEN p)
872 96320 : { return gen_factorback(L, e, &Fpmul, &Fppow, (void*)p); }
873 :
874 : GEN
875 15665411 : factorback2(GEN L, GEN e) { return gen_factorback(L, e, &mul, &powi, NULL); }
876 : GEN
877 1180438 : factorback(GEN fa) { return factorback2(fa, NULL); }
878 :
879 : GEN
880 14 : vecprod(GEN v)
881 : {
882 14 : pari_sp av = avma;
883 14 : if (!is_vec_t(typ(v)))
884 0 : pari_err_TYPE("vecprod", v);
885 14 : if (lg(v) == 1) return gen_1;
886 7 : return gerepilecopy(av, gen_product(v, NULL, mul));
887 : }
888 :
889 : static int
890 63 : RgX_is_irred_i(GEN x)
891 : {
892 : GEN y, p, pol;
893 63 : long l = lg(x), pa;
894 :
895 63 : if (!signe(x) || l <= 3) return 0;
896 63 : switch(RgX_type(x,&p,&pol,&pa))
897 : {
898 14 : case t_INTMOD: return FpX_is_irred(RgX_to_FpX(x,p), p);
899 0 : case t_COMPLEX: return l == 4;
900 : case t_REAL:
901 0 : if (l == 4) return 1;
902 0 : if (l > 5) return 0;
903 0 : return gsigne(RgX_disc(x)) > 0;
904 : }
905 49 : y = factor(x);
906 49 : return (lg(gcoeff(y,1,1))==l);
907 : }
908 : static int
909 63 : RgX_is_irred(GEN x)
910 : {
911 63 : pari_sp av = avma;
912 63 : int r = RgX_is_irred_i(x);
913 63 : avma = av; return r;
914 : }
915 : long
916 63 : isirreducible(GEN x)
917 : {
918 63 : switch(typ(x))
919 : {
920 0 : case t_INT: case t_REAL: case t_FRAC: return 0;
921 63 : case t_POL: return RgX_is_irred(x);
922 : }
923 0 : pari_err_TYPE("isirreducible",x);
924 0 : return 0;
925 : }
926 :
927 : /*******************************************************************/
928 : /* */
929 : /* GENERIC GCD */
930 : /* */
931 : /*******************************************************************/
932 : /* x is a COMPLEX or a QUAD */
933 : static GEN
934 357 : triv_cont_gcd(GEN x, GEN y)
935 : {
936 357 : pari_sp av = avma;
937 : GEN c;
938 357 : if (typ(x)==t_COMPLEX)
939 : {
940 35 : GEN a = gel(x,1), b = gel(x,2);
941 35 : if (typ(a) == t_REAL || typ(b) == t_REAL) return gen_1;
942 21 : c = ggcd(a,b);
943 : }
944 : else
945 322 : c = ggcd(gel(x,2),gel(x,3));
946 343 : return gerepileupto(av, ggcd(c,y));
947 : }
948 :
949 : /* y is a PADIC, x a rational number or an INTMOD */
950 : static GEN
951 182 : padic_gcd(GEN x, GEN y)
952 : {
953 182 : GEN p = gel(y,2);
954 182 : long v = gvaluation(x,p), w = valp(y);
955 182 : if (w < v) v = w;
956 182 : return powis(p, v);
957 : }
958 :
959 : /* x,y in Z[i], at least one of which is t_COMPLEX */
960 : static GEN
961 28 : gauss_gcd(GEN x, GEN y)
962 : {
963 28 : pari_sp av = avma;
964 : GEN dx, dy;
965 28 : dx = denom(x); x = gmul(x, dx);
966 28 : dy = denom(y); y = gmul(y, dy);
967 77 : while (!gequal0(y))
968 : {
969 21 : GEN z = gsub(x, gmul(ground(gdiv(x,y)), y));
970 21 : x = y; y = z;
971 : }
972 28 : x = gauss_normal(x);
973 28 : if (typ(x) == t_COMPLEX)
974 : {
975 21 : if (gequal0(gel(x,2))) x = gel(x,1);
976 21 : else if (gequal0(gel(x,1))) x = gel(x,2);
977 : }
978 28 : return gerepileupto(av, gdiv(x, lcmii(dx, dy)));
979 : }
980 :
981 : static int
982 35 : c_is_rational(GEN x)
983 35 : { return is_rational_t(typ(gel(x,1))) && is_rational_t(typ(gel(x,2))); }
984 : static GEN
985 49 : c_zero_gcd(GEN c)
986 : {
987 49 : GEN x = gel(c,1), y = gel(c,2);
988 49 : long tx = typ(x), ty = typ(y);
989 49 : if (tx == t_REAL || ty == t_REAL) return gen_1;
990 21 : if (tx == t_PADIC || tx == t_INTMOD
991 14 : || ty == t_PADIC || ty == t_INTMOD) return ggcd(x, y);
992 14 : return gauss_gcd(c, gen_0);
993 : }
994 :
995 : /* gcd(x, 0) */
996 : static GEN
997 8192154 : zero_gcd(GEN x)
998 : {
999 : pari_sp av;
1000 8192154 : switch(typ(x))
1001 : {
1002 35528 : case t_INT: return absi(x);
1003 1344 : case t_FRAC: return absfrac(x);
1004 49 : case t_COMPLEX: return c_zero_gcd(x);
1005 2016 : case t_REAL: return gen_1;
1006 56 : case t_PADIC: return powis(gel(x,2), valp(x));
1007 210 : case t_SER: return pol_xn(valp(x), varn(x));
1008 : case t_POLMOD: {
1009 13966 : GEN d = gel(x,2);
1010 13966 : if (typ(d) == t_POL && varn(d) == varn(gel(x,1))) return content(d);
1011 497 : return isinexact(d)? zero_gcd(d): gcopy(d);
1012 : }
1013 : case t_POL:
1014 7921813 : if (!isinexact(x)) break;
1015 0 : av = avma;
1016 0 : return gerepileupto(av,
1017 0 : monomialcopy(content(x), RgX_val(x), varn(x))
1018 : );
1019 :
1020 : case t_RFRAC:
1021 216969 : if (!isinexact(x)) break;
1022 0 : av = avma;
1023 0 : return gerepileupto(av,
1024 0 : gdiv(zero_gcd(gel(x,1)), gel(x,2))
1025 : );
1026 : }
1027 8138985 : return gcopy(x);
1028 : }
1029 : /* z is an exact zero, t_INT, t_INTMOD or t_FFELT */
1030 : static GEN
1031 8177022 : zero_gcd2(GEN y, GEN z)
1032 : {
1033 : pari_sp av;
1034 8177022 : switch(typ(z))
1035 : {
1036 8172745 : case t_INT: return zero_gcd(y);
1037 : case t_INTMOD:
1038 3367 : av = avma;
1039 3367 : return gerepileupto(av, gmul(y, mkintmod(gen_1,gel(z,1))));
1040 : case t_FFELT:
1041 910 : av = avma;
1042 910 : return gerepileupto(av, gmul(y, FF_1(z)));
1043 : default:
1044 0 : pari_err_TYPE("zero_gcd", z);
1045 : }
1046 0 : return NULL;
1047 : }
1048 : /* tx = t_RFRAC, y considered as constant */
1049 : static GEN
1050 864892 : cont_gcd_rfrac(GEN x, GEN y)
1051 : {
1052 864892 : pari_sp av = avma;
1053 864892 : GEN cx; x = primitive_part(x, &cx);
1054 864892 : return gerepileupto(av, gred_rfrac_simple(ggcd(cx? cx: gen_1, y), gel(x,2)));
1055 : }
1056 : /* tx = t_POL, y considered as constant */
1057 : static GEN
1058 2271783 : cont_gcd_pol(GEN x, GEN y)
1059 : {
1060 2271783 : pari_sp av = avma;
1061 2271783 : return gerepileupto(av, scalarpol(ggcd(content(x),y), varn(x)));
1062 : }
1063 : /* !is_const_t(tx), tx != t_POL,t_RFRAC, y considered as constant */
1064 : static GEN
1065 13322 : cont_gcd_gen(GEN x, GEN y)
1066 : {
1067 13322 : pari_sp av = avma;
1068 13322 : return gerepileupto(av, ggcd(content(x),y));
1069 : }
1070 : /* !is_const(tx), y considered as constant */
1071 : static GEN
1072 3093494 : cont_gcd(GEN x, long tx, GEN y)
1073 : {
1074 3093494 : switch(tx)
1075 : {
1076 808445 : case t_RFRAC: return cont_gcd_rfrac(x,y);
1077 2271727 : case t_POL: return cont_gcd_pol(x,y);
1078 13322 : default: return cont_gcd_gen(x,y);
1079 : }
1080 : }
1081 : static GEN
1082 699802 : gcdiq(GEN x, GEN y)
1083 : {
1084 : GEN z;
1085 699802 : if (!signe(x)) return Q_abs(y);
1086 212462 : z = cgetg(3,t_FRAC);
1087 212462 : gel(z,1) = gcdii(x,gel(y,1));
1088 212462 : gel(z,2) = icopy(gel(y,2));
1089 212462 : return z;
1090 : }
1091 : static GEN
1092 673892 : gcdqq(GEN x, GEN y)
1093 : {
1094 673892 : GEN z = cgetg(3,t_FRAC);
1095 673892 : gel(z,1) = gcdii(gel(x,1), gel(y,1));
1096 673892 : gel(z,2) = lcmii(gel(x,2), gel(y,2));
1097 673892 : return z;
1098 : }
1099 : /* assume x,y t_INT or t_FRAC */
1100 : GEN
1101 51421232 : Q_gcd(GEN x, GEN y)
1102 : {
1103 51421232 : long tx = typ(x), ty = typ(y);
1104 51421232 : if (tx == t_INT)
1105 50206759 : { return (ty == t_INT)? gcdii(x,y): gcdiq(x,y); }
1106 : else
1107 1214473 : { return (ty == t_INT)? gcdiq(y,x): gcdqq(x,y); }
1108 : }
1109 :
1110 : GEN
1111 21469559 : ggcd(GEN x, GEN y)
1112 : {
1113 21469559 : long vx, vy, tx = typ(x), ty = typ(y);
1114 : pari_sp av, tetpil;
1115 : GEN p1,z;
1116 :
1117 42939118 : if (is_noncalc_t(tx) || is_matvec_t(tx) ||
1118 42939118 : is_noncalc_t(ty) || is_matvec_t(ty)) pari_err_TYPE2("gcd",x,y);
1119 21469559 : if (tx>ty) { swap(x,y); lswap(tx,ty); }
1120 : /* tx <= ty */
1121 21469559 : z = gisexactzero(x); if (z) return zero_gcd2(y,z);
1122 18447013 : z = gisexactzero(y); if (z) return zero_gcd2(x,z);
1123 13292537 : if (is_const_t(tx))
1124 : {
1125 6365476 : if (ty == tx) switch(tx)
1126 : {
1127 : case t_INT:
1128 3245056 : return gcdii(x,y);
1129 :
1130 1211 : case t_INTMOD: z=cgetg(3,t_INTMOD);
1131 1211 : if (equalii(gel(x,1),gel(y,1)))
1132 1211 : gel(z,1) = icopy(gel(x,1));
1133 : else
1134 0 : gel(z,1) = gcdii(gel(x,1),gel(y,1));
1135 1211 : if (gequal1(gel(z,1))) gel(z,2) = gen_0;
1136 : else
1137 : {
1138 1211 : av = avma; p1 = gcdii(gel(z,1),gel(x,2));
1139 1211 : if (!is_pm1(p1))
1140 : {
1141 0 : p1 = gcdii(p1,gel(y,2));
1142 0 : if (equalii(p1, gel(z,1))) { cgiv(p1); p1 = gen_0; }
1143 0 : else p1 = gerepileuptoint(av, p1);
1144 : }
1145 1211 : gel(z,2) = p1;
1146 : }
1147 1211 : return z;
1148 :
1149 : case t_FRAC:
1150 34785 : return gcdqq(x,y);
1151 :
1152 : case t_FFELT:
1153 7931 : if (!FF_samefield(x,y)) pari_err_OP("gcd",x,y);
1154 7931 : return FF_equal0(x) && FF_equal0(y)? FF_zero(y): FF_1(y);
1155 :
1156 : case t_COMPLEX:
1157 14 : if (c_is_rational(x) && c_is_rational(y)) return gauss_gcd(x,y);
1158 7 : return triv_cont_gcd(y,x);
1159 :
1160 : case t_PADIC:
1161 14 : if (!equalii(gel(x,2),gel(y,2))) return gen_1;
1162 7 : return powis(gel(y,2), minss(valp(x), valp(y)));
1163 :
1164 : case t_QUAD:
1165 133 : av=avma; p1=gdiv(x,y);
1166 133 : if (gequal0(gel(p1,3)))
1167 : {
1168 14 : p1=gel(p1,2);
1169 14 : if (typ(p1)==t_INT) { avma=av; return gcopy(y); }
1170 7 : tetpil=avma; return gerepile(av,tetpil, gdiv(y,gel(p1,2)));
1171 : }
1172 119 : if (typ(gel(p1,2))==t_INT && typ(gel(p1,3))==t_INT) {avma=av; return gcopy(y);}
1173 112 : p1 = ginv(p1); avma=av;
1174 112 : if (typ(gel(p1,2))==t_INT && typ(gel(p1,3))==t_INT) return gcopy(x);
1175 105 : return triv_cont_gcd(y,x);
1176 :
1177 0 : default: return gen_1; /* t_REAL */
1178 : }
1179 3076332 : if (is_const_t(ty)) switch(tx)
1180 : {
1181 : case t_INT:
1182 10754 : switch(ty)
1183 : {
1184 42 : case t_INTMOD: z = cgetg(3,t_INTMOD);
1185 42 : gel(z,1) = icopy(gel(y,1)); av = avma;
1186 42 : p1 = gcdii(gel(y,1),gel(y,2));
1187 42 : if (!is_pm1(p1)) {
1188 14 : p1 = gcdii(x,p1);
1189 14 : if (equalii(p1, gel(z,1))) { cgiv(p1); p1 = gen_0; }
1190 : else
1191 14 : p1 = gerepileuptoint(av, p1);
1192 : }
1193 42 : gel(z,2) = p1; return z;
1194 :
1195 5929 : case t_REAL: return gen_1;
1196 :
1197 : case t_FRAC:
1198 3362 : return gcdiq(x,y);
1199 :
1200 : case t_COMPLEX:
1201 0 : if (c_is_rational(y)) return gauss_gcd(x,y);
1202 0 : return triv_cont_gcd(y,x);
1203 :
1204 : case t_FFELT:
1205 1127 : if (!FF_equal0(y)) return FF_1(y);
1206 0 : return dvdii(x, gel(y,4))? FF_zero(y): FF_1(y);
1207 :
1208 : case t_PADIC:
1209 168 : return padic_gcd(x,y);
1210 :
1211 : case t_QUAD:
1212 126 : return triv_cont_gcd(y,x);
1213 : default:
1214 0 : pari_err_TYPE2("gcd",x,y);
1215 : }
1216 :
1217 : case t_REAL:
1218 14 : switch(ty)
1219 : {
1220 : case t_INTMOD:
1221 : case t_FFELT:
1222 14 : case t_PADIC: pari_err_TYPE2("gcd",x,y);
1223 0 : default: return gen_1;
1224 : }
1225 :
1226 : case t_INTMOD:
1227 63 : switch(ty)
1228 : {
1229 : case t_FRAC:
1230 14 : av = avma; p1=gcdii(gel(x,1),gel(y,2)); avma = av;
1231 14 : if (!is_pm1(p1)) pari_err_OP("gcd",x,y);
1232 7 : return ggcd(gel(y,1), x);
1233 :
1234 : case t_FFELT:
1235 : {
1236 28 : GEN p = gel(y,4);
1237 28 : if (!dvdii(gel(x,1), p)) pari_err_OP("gcd",x,y);
1238 21 : if (!FF_equal0(y)) return FF_1(y);
1239 0 : return dvdii(gel(x,2),p)? FF_zero(y): FF_1(y);
1240 : }
1241 :
1242 : case t_COMPLEX: case t_QUAD:
1243 14 : return triv_cont_gcd(y,x);
1244 :
1245 : case t_PADIC:
1246 7 : return padic_gcd(x,y);
1247 :
1248 0 : default: pari_err_TYPE2("gcd",x,y);
1249 : }
1250 :
1251 : case t_FRAC:
1252 140 : switch(ty)
1253 : {
1254 : case t_COMPLEX:
1255 14 : if (c_is_rational(y)) return gauss_gcd(x,y);
1256 : case t_QUAD:
1257 84 : return triv_cont_gcd(y,x);
1258 : case t_FFELT:
1259 : {
1260 42 : GEN p = gel(y,4);
1261 42 : if (dvdii(gel(x,2), p)) pari_err_OP("gcd",x,y);
1262 21 : if (!FF_equal0(y)) return FF_1(y);
1263 0 : return dvdii(gel(x,1),p)? FF_zero(y): FF_1(y);
1264 : }
1265 :
1266 : case t_PADIC:
1267 7 : return padic_gcd(x,y);
1268 :
1269 0 : default: pari_err_TYPE2("gcd",x,y);
1270 : }
1271 : case t_FFELT:
1272 70 : switch(ty)
1273 : {
1274 : case t_PADIC:
1275 : {
1276 42 : GEN p = gel(y,2);
1277 42 : long v = valp(y);
1278 42 : if (!equalii(p, gel(x,4)) || v < 0) pari_err_OP("gcd",x,y);
1279 14 : return (v && FF_equal0(x))? FF_zero(x): FF_1(x);
1280 : }
1281 28 : default: pari_err_TYPE2("gcd",x,y);
1282 : }
1283 :
1284 : case t_COMPLEX:
1285 14 : switch(ty)
1286 : {
1287 : case t_PADIC:
1288 14 : case t_QUAD: return triv_cont_gcd(x,y);
1289 0 : default: pari_err_TYPE2("gcd",x,y);
1290 : }
1291 :
1292 : case t_PADIC:
1293 7 : switch(ty)
1294 : {
1295 7 : case t_QUAD: return triv_cont_gcd(y,x);
1296 0 : default: pari_err_TYPE2("gcd",x,y);
1297 : }
1298 :
1299 0 : default: return gen_1; /* tx = t_REAL */
1300 : }
1301 3065270 : return cont_gcd(y,ty, x);
1302 : }
1303 :
1304 6927061 : if (tx == t_POLMOD)
1305 : {
1306 13394 : if (ty == t_POLMOD)
1307 : {
1308 13051 : GEN T = gel(x,1);
1309 13051 : z = cgetg(3,t_POLMOD);
1310 13051 : T = RgX_equal_var(T,gel(y,1))? RgX_copy(T): RgX_gcd(T, gel(y,1));
1311 13051 : gel(z,1) = T;
1312 13051 : if (degpol(T) <= 0) gel(z,2) = gen_0;
1313 : else
1314 : {
1315 : GEN X, Y, d;
1316 13051 : av = avma; X = gel(x,2); Y = gel(y,2);
1317 13051 : d = ggcd(content(X), content(Y));
1318 13051 : if (!gequal1(d)) { X = gdiv(X,d); Y = gdiv(Y,d); }
1319 13051 : p1 = ggcd(T, X);
1320 13051 : gel(z,2) = gerepileupto(av, gmul(d, ggcd(p1, Y)));
1321 : }
1322 13051 : return z;
1323 : }
1324 343 : vx = varn(gel(x,1));
1325 343 : switch(ty)
1326 : {
1327 : case t_POL:
1328 315 : vy = varn(y);
1329 315 : if (varncmp(vy,vx) < 0) return cont_gcd_pol(y, x);
1330 259 : z = cgetg(3,t_POLMOD);
1331 259 : gel(z,1) = RgX_copy(gel(x,1));
1332 259 : av = avma; p1 = ggcd(gel(x,1),gel(x,2));
1333 259 : gel(z,2) = gerepileupto(av, ggcd(p1,y));
1334 259 : return z;
1335 :
1336 : case t_RFRAC:
1337 28 : vy = varn(gel(y,2));
1338 28 : if (varncmp(vy,vx) < 0) return cont_gcd_rfrac(y, x);
1339 28 : av = avma;
1340 28 : p1 = ggcd(gel(x,1),gel(y,2));
1341 28 : if (degpol(p1)) pari_err_OP("gcd",x,y);
1342 21 : avma = av; return gdiv(ggcd(gel(y,1),x), content(gel(y,2)));
1343 : }
1344 : }
1345 :
1346 6913667 : vx = gvar(x);
1347 6913667 : vy = gvar(y);
1348 6913667 : if (varncmp(vy, vx) < 0) return cont_gcd(y,ty, x);
1349 6901571 : if (varncmp(vy, vx) > 0) return cont_gcd(x,tx, y);
1350 :
1351 : /* vx = vy: same main variable */
1352 6885443 : switch(tx)
1353 : {
1354 : case t_POL:
1355 6724537 : switch(ty)
1356 : {
1357 6668076 : case t_POL: return RgX_gcd(x,y);
1358 : case t_SER:
1359 14 : z = ggcd(content(x), content(y));
1360 14 : return monomialcopy(z, minss(valp(y),gval(x,vx)), vx);
1361 56447 : case t_RFRAC: return cont_gcd_rfrac(y, x);
1362 : }
1363 0 : break;
1364 :
1365 : case t_SER:
1366 14 : z = ggcd(content(x), content(y));
1367 14 : switch(ty)
1368 : {
1369 7 : case t_SER: return monomialcopy(z, minss(valp(x),valp(y)), vx);
1370 7 : case t_RFRAC: return monomialcopy(z, minss(valp(x),gval(y,vx)), vx);
1371 : }
1372 0 : break;
1373 :
1374 : case t_RFRAC:
1375 : {
1376 160892 : GEN xd = gel(x,2), yd = gel(y,2);
1377 160892 : if (ty != t_RFRAC) pari_err_TYPE2("gcd",x,y);
1378 160892 : z = cgetg(3,t_RFRAC); av = avma;
1379 160892 : gel(z,2) = gerepileupto(av, RgX_mul(xd, RgX_div(yd, RgX_gcd(xd, yd))));
1380 160892 : gel(z,1) = ggcd(gel(x,1), gel(y,1)); return z;
1381 : }
1382 : }
1383 0 : pari_err_TYPE2("gcd",x,y);
1384 : return NULL; /* LCOV_EXCL_LINE */
1385 : }
1386 : GEN
1387 3960 : ggcd0(GEN x, GEN y) { return y? ggcd(x,y): content(x); }
1388 :
1389 : static GEN
1390 133 : fix_lcm(GEN x)
1391 : {
1392 : GEN t;
1393 133 : switch(typ(x))
1394 : {
1395 35 : case t_INT: if (signe(x)<0) x = negi(x);
1396 35 : break;
1397 : case t_POL:
1398 98 : if (lg(x) <= 2) break;
1399 98 : t = leading_coeff(x);
1400 98 : if (typ(t) == t_INT && signe(t) < 0) x = gneg(x);
1401 : }
1402 133 : return x;
1403 : }
1404 : GEN
1405 2891 : glcm0(GEN x, GEN y)
1406 : {
1407 2891 : if (!y) return fix_lcm(gassoc_proto(glcm,x,y));
1408 2807 : return glcm(x,y);
1409 : }
1410 : GEN
1411 3276 : glcm(GEN x, GEN y)
1412 : {
1413 : pari_sp av;
1414 : GEN z;
1415 3276 : if (typ(x)==t_INT && typ(y)==t_INT) return lcmii(x,y);
1416 49 : av = avma;
1417 49 : z = ggcd(x,y); if (!gequal1(z)) y = gdiv(y,z);
1418 49 : return gerepileupto(av, fix_lcm(gmul(x,y)));
1419 : }
1420 :
1421 : /* x + r ~ x ? Assume x,r are t_POL, deg(r) <= deg(x) */
1422 : static int
1423 2107 : pol_approx0(GEN r, GEN x, int exact)
1424 : {
1425 : long i, lx,lr;
1426 2107 : if (exact) return gequal0(r);
1427 0 : lx = lg(x);
1428 0 : lr = lg(r); if (lr < lx) lx = lr;
1429 0 : for (i=2; i<lx; i++)
1430 0 : if (!approx_0(gel(r,i), gel(x,i))) return 0;
1431 0 : return 1;
1432 : }
1433 :
1434 : GEN
1435 777 : RgX_gcd_simple(GEN x, GEN y)
1436 : {
1437 777 : pari_sp av1, av = avma;
1438 777 : GEN r, yorig = y;
1439 777 : int exact = !(isinexactreal(x) || isinexactreal(y));
1440 :
1441 : for(;;)
1442 : {
1443 2107 : av1 = avma; r = RgX_rem(x,y);
1444 2107 : if (pol_approx0(r, x, exact))
1445 : {
1446 777 : avma = av1;
1447 777 : if (y == yorig) return RgX_copy(y);
1448 665 : y = normalizepol_approx(y, lg(y));
1449 665 : if (lg(y) == 3) { avma = av; return pol_1(varn(x)); }
1450 98 : return gerepileupto(av,y);
1451 : }
1452 1330 : x = y; y = r;
1453 1330 : if (gc_needed(av,1)) {
1454 0 : if(DEBUGMEM>1) pari_warn(warnmem,"RgX_gcd_simple");
1455 0 : gerepileall(av,2, &x,&y);
1456 : }
1457 1330 : }
1458 : }
1459 : GEN
1460 0 : RgX_extgcd_simple(GEN a, GEN b, GEN *pu, GEN *pv)
1461 : {
1462 0 : pari_sp av = avma;
1463 : GEN q, r, d, d1, u, v, v1;
1464 0 : int exact = !(isinexactreal(a) || isinexactreal(b));
1465 :
1466 0 : d = a; d1 = b; v = gen_0; v1 = gen_1;
1467 : for(;;)
1468 : {
1469 0 : if (pol_approx0(d1, a, exact)) break;
1470 0 : q = poldivrem(d,d1, &r);
1471 0 : v = gsub(v, gmul(q,v1));
1472 0 : u=v; v=v1; v1=u;
1473 0 : u=r; d=d1; d1=u;
1474 0 : }
1475 0 : u = gsub(d, gmul(b,v));
1476 0 : u = RgX_div(u,a);
1477 :
1478 0 : gerepileall(av, 3, &u,&v,&d);
1479 0 : *pu = u;
1480 0 : *pv = v; return d;
1481 : }
1482 : /*******************************************************************/
1483 : /* */
1484 : /* CONTENT / PRIMITIVE PART */
1485 : /* */
1486 : /*******************************************************************/
1487 :
1488 : GEN
1489 77075147 : content(GEN x)
1490 : {
1491 77075147 : long lx, i, t, tx = typ(x);
1492 77075147 : pari_sp av = avma;
1493 : GEN c;
1494 :
1495 77075147 : if (is_scalar_t(tx)) return zero_gcd(x);
1496 77058027 : switch(tx)
1497 : {
1498 : case t_RFRAC:
1499 : {
1500 864906 : GEN n = gel(x,1), d = gel(x,2);
1501 : /* -- varncmp(vn, vd) < 0 can't happen
1502 : * -- if n is POLMOD, its main variable (in the sense of gvar2)
1503 : * has lower priority than denominator */
1504 864906 : if (typ(n) == t_POLMOD || varncmp(gvar(n), varn(d)) > 0)
1505 824650 : n = isinexact(n)? zero_gcd(n): gcopy(n);
1506 : else
1507 40256 : n = content(n);
1508 864906 : return gerepileupto(av, gdiv(n, content(d)));
1509 : }
1510 :
1511 : case t_VEC: case t_COL:
1512 7546024 : lx = lg(x); if (lx==1) return gen_0;
1513 7546017 : break;
1514 :
1515 : case t_MAT:
1516 : {
1517 : long hx, j;
1518 238 : lx = lg(x);
1519 238 : if (lx == 1) return gen_0;
1520 231 : hx = lgcols(x);
1521 231 : if (hx == 1) return gen_0;
1522 224 : if (lx == 2) { x = gel(x,1); lx = lg(x); break; }
1523 224 : if (hx == 2) { x = row_i(x, 1, 1, lx-1); break; }
1524 224 : c = content(gel(x,1));
1525 448 : for (j=2; j<lx; j++)
1526 224 : for (i=1; i<hx; i++) c = ggcd(c,gcoeff(x,i,j));
1527 224 : if (typ(c) == t_INTMOD || isinexact(c)) { avma=av; return gen_1; }
1528 224 : return gerepileupto(av,c);
1529 : }
1530 :
1531 : case t_POL: case t_SER:
1532 68646824 : lx = lg(x); if (lx == 2) return gen_0;
1533 68646810 : break;
1534 21 : case t_VECSMALL: return utoi(zv_content(x));
1535 : case t_QFR: case t_QFI:
1536 14 : lx = 4; break;
1537 :
1538 0 : default: pari_err_TYPE("content",x);
1539 : return NULL; /* LCOV_EXCL_LINE */
1540 : }
1541 224883062 : for (i=lontyp[tx]; i<lx; i++)
1542 161336515 : if (typ(gel(x,i)) != t_INT) break;
1543 76192841 : lx--; c = gel(x,lx);
1544 76192841 : t = typ(c); if (is_matvec_t(t)) c = content(c);
1545 76192841 : if (i > lx)
1546 : { /* integer coeffs */
1547 128310843 : while (lx-- > lontyp[tx])
1548 : {
1549 64015524 : c = gcdii(c, gel(x,lx));
1550 64015524 : if (is_pm1(c)) { avma=av; return gen_1; }
1551 : }
1552 : }
1553 : else
1554 : {
1555 12646294 : if (isinexact(c)) c = zero_gcd(c);
1556 43402706 : while (lx-- > lontyp[tx])
1557 : {
1558 18110118 : GEN d = gel(x,lx);
1559 18110118 : t = typ(d); if (is_matvec_t(t)) d = content(d);
1560 18110118 : c = ggcd(c, d);
1561 : }
1562 12646294 : if (isinexact(c)) { avma=av; return gen_1; }
1563 : }
1564 13395066 : switch(typ(c))
1565 : {
1566 : case t_INT:
1567 757627 : if (signe(c) < 0) c = negi(c);
1568 757627 : break;
1569 : case t_VEC: case t_COL: case t_MAT:
1570 0 : pari_err_TYPE("content",x);
1571 : }
1572 :
1573 13395066 : return av==avma? gcopy(c): gerepileupto(av,c);
1574 : }
1575 :
1576 : GEN
1577 1085271 : primitive_part(GEN x, GEN *ptc)
1578 : {
1579 1085271 : pari_sp av = avma;
1580 1085271 : GEN c = content(x);
1581 1085271 : if (gequal1(c)) { avma = av; c = NULL; }
1582 35546 : else if (!gequal0(c)) x = gdiv(x,c);
1583 1085271 : if (ptc) *ptc = c;
1584 1085271 : return x;
1585 : }
1586 : GEN
1587 2485 : primpart(GEN x) { return primitive_part(x, NULL); }
1588 :
1589 : /* As content(), but over Q. Treats polynomial as elts of Q[x1,...xn], instead
1590 : * of Q(x2,...,xn)[x1] */
1591 : GEN
1592 73637092 : Q_content(GEN x)
1593 : {
1594 : long i, l;
1595 : GEN d;
1596 : pari_sp av;
1597 :
1598 73637092 : switch(typ(x))
1599 : {
1600 57509060 : case t_INT: return absi(x);
1601 858124 : case t_FRAC: return absfrac(x);
1602 :
1603 : case t_VEC: case t_COL: case t_MAT:
1604 7452536 : l = lg(x); if (l == 1) return gen_1;
1605 7452487 : av = avma; d = Q_content(gel(x,1));
1606 33831456 : for (i=2; i<l; i++)
1607 : {
1608 26378969 : d = Q_gcd(d, Q_content(gel(x,i)));
1609 26378969 : if ((i & 255) == 0) d = gerepileupto(av, d);
1610 : }
1611 7452487 : return gerepileupto(av, d);
1612 :
1613 : case t_POL:
1614 7812262 : l = lg(x); if (l == 2) return gen_0;
1615 7722732 : av = avma; d = Q_content(gel(x,2));
1616 7722732 : for (i=3; i<l; i++) d = Q_gcd(d, Q_content(gel(x,i)));
1617 7722732 : return gerepileupto(av, d);
1618 5089 : case t_POLMOD: return Q_content(gel(x,2));
1619 21 : case t_COMPLEX: return Q_gcd(Q_content(gel(x,1)), Q_content(gel(x,2)));
1620 : }
1621 0 : pari_err_TYPE("Q_content",x);
1622 : return NULL; /* LCOV_EXCL_LINE */
1623 : }
1624 :
1625 : GEN
1626 13153 : ZX_content(GEN x)
1627 : {
1628 13153 : long i, l = lg(x);
1629 : GEN d;
1630 : pari_sp av;
1631 :
1632 13153 : if (l == 2) return gen_0;
1633 13153 : d = gel(x,2);
1634 13153 : if (l == 3) return absi(d);
1635 9177 : av = avma;
1636 9177 : for (i=3; !is_pm1(d) && i<l; i++) d = gcdii(d, gel(x,i));
1637 9177 : if (signe(d) < 0) d = absi(d);
1638 9177 : return gerepileuptoint(av, d);
1639 : }
1640 :
1641 : /* NOT MEMORY CLEAN (because of t_FRAC).
1642 : * As denom(), but over Q. Treats polynomial as elts of Q[x1,...xn], instead
1643 : * of Q(x2,...,xn)[x1] */
1644 : GEN
1645 56458602 : Q_denom(GEN x)
1646 : {
1647 : long i, l;
1648 : GEN d, D;
1649 : pari_sp av;
1650 :
1651 56458602 : switch(typ(x))
1652 : {
1653 35365282 : case t_INT: return gen_1;
1654 10138768 : case t_FRAC: return gel(x,2);
1655 :
1656 : case t_VEC: case t_COL: case t_MAT:
1657 8197906 : l = lg(x); if (l == 1) return gen_1;
1658 8132155 : av = avma; d = Q_denom(gel(x,1));
1659 35847881 : for (i=2; i<l; i++)
1660 : {
1661 27715726 : D = Q_denom(gel(x,i));
1662 27715726 : if (D != gen_1) d = lcmii(d, D);
1663 27715726 : if ((i & 255) == 0) d = gerepileuptoint(av, d);
1664 : }
1665 8132155 : return gerepileuptoint(av, d);
1666 :
1667 : case t_POL:
1668 2705371 : l = lg(x); if (l == 2) return gen_1;
1669 2659500 : av = avma; d = Q_denom(gel(x,2));
1670 7691113 : for (i=3; i<l; i++)
1671 : {
1672 5031634 : D = Q_denom(gel(x,i));
1673 5031634 : if (D != gen_1) d = lcmii(d, D);
1674 : }
1675 2659479 : return gerepileuptoint(av, d);
1676 :
1677 51254 : case t_POLMOD: return Q_denom(gel(x,2));
1678 : }
1679 21 : pari_err_TYPE("Q_denom",x);
1680 : return NULL; /* LCOV_EXCL_LINE */
1681 : }
1682 :
1683 : GEN
1684 11296720 : Q_remove_denom(GEN x, GEN *ptd)
1685 : {
1686 11296720 : GEN d = Q_denom(x);
1687 11296699 : if (d == gen_1) d = NULL; else x = Q_muli_to_int(x,d);
1688 11296699 : if (ptd) *ptd = d;
1689 11296699 : return x;
1690 : }
1691 :
1692 : /* return y = x * d, assuming x rational, and d,y integral */
1693 : GEN
1694 21027217 : Q_muli_to_int(GEN x, GEN d)
1695 : {
1696 : long i, l;
1697 : GEN y, xn, xd;
1698 : pari_sp av;
1699 :
1700 21027217 : if (typ(d) != t_INT) pari_err_TYPE("Q_muli_to_int",d);
1701 21027217 : switch (typ(x))
1702 : {
1703 : case t_INT:
1704 6914785 : return mulii(x,d);
1705 :
1706 : case t_FRAC:
1707 9667687 : xn = gel(x,1);
1708 9667687 : xd = gel(x,2); av = avma;
1709 9667687 : y = mulii(xn, diviiexact(d, xd));
1710 9667687 : return gerepileuptoint(av, y);
1711 :
1712 : case t_VEC: case t_COL: case t_MAT:
1713 3884294 : y = cgetg_copy(x, &l);
1714 3884294 : for (i=1; i<l; i++) gel(y,i) = Q_muli_to_int(gel(x,i), d);
1715 3884294 : return y;
1716 :
1717 : case t_POL:
1718 544463 : y = cgetg_copy(x, &l); y[1] = x[1];
1719 544463 : for (i=2; i<l; i++) gel(y,i) = Q_muli_to_int(gel(x,i), d);
1720 544463 : return y;
1721 :
1722 : case t_POLMOD:
1723 15988 : y = cgetg(3, t_POLMOD);
1724 15988 : gel(y,1) = RgX_copy(gel(x,1));
1725 15988 : gel(y,2) = Q_muli_to_int(gel(x,2), d);
1726 15988 : return y;
1727 : }
1728 0 : pari_err_TYPE("Q_muli_to_int",x);
1729 : return NULL; /* LCOV_EXCL_LINE */
1730 : }
1731 :
1732 : static void
1733 1818283 : rescale_init(GEN c, int *exact, long *emin, GEN *D)
1734 : {
1735 : long e;
1736 1818283 : switch(typ(c))
1737 : {
1738 : case t_REAL:
1739 1216376 : *exact = 0;
1740 1216376 : if (!signe(c)) return;
1741 1214126 : e = expo(c) - bit_prec(c);
1742 1214126 : break;
1743 : case t_INT:
1744 601809 : if (!signe(c)) return;
1745 247534 : e = expi(c) + 32;
1746 247534 : break;
1747 : case t_FRAC:
1748 70 : e = expi(gel(c,1)) - expi(gel(c,2)) + 32;
1749 70 : if (exact) *D = lcmii(*D, gel(c,2));
1750 70 : break;
1751 : default:
1752 28 : pari_err_TYPE("rescale_to_int",c);
1753 : return; /* LCOV_EXCL_LINE */
1754 : }
1755 1461730 : if (e < *emin) *emin = e;
1756 : }
1757 : GEN
1758 222557 : RgM_rescale_to_int(GEN x)
1759 : {
1760 222557 : long lx = lg(x), i,j, hx, emin;
1761 : GEN D;
1762 : int exact;
1763 :
1764 222557 : if (lx == 1) return cgetg(1,t_MAT);
1765 222557 : hx = lgcols(x);
1766 222557 : exact = 1;
1767 222557 : emin = HIGHEXPOBIT;
1768 222557 : D = gen_1;
1769 805816 : for (j = 1; j < lx; j++)
1770 583287 : for (i = 1; i < hx; i++) rescale_init(gcoeff(x,i,j), &exact, &emin, &D);
1771 222529 : if (exact) return D == gen_1 ? x: Q_muli_to_int(x, D);
1772 222473 : return grndtoi(gmul2n(x, -emin), &i);
1773 : }
1774 : GEN
1775 238 : RgX_rescale_to_int(GEN x)
1776 : {
1777 238 : long lx = lg(x), i, emin;
1778 : GEN D;
1779 : int exact;
1780 238 : if (lx == 2) return gcopy(x); /* rare */
1781 238 : exact = 1;
1782 238 : emin = HIGHEXPOBIT;
1783 238 : D = gen_1;
1784 238 : for (i = 2; i < lx; i++) rescale_init(gel(x,i), &exact, &emin, &D);
1785 238 : if (exact) return D == gen_1 ? x: Q_muli_to_int(x, D);
1786 238 : return grndtoi(gmul2n(x, -emin), &i);
1787 : }
1788 :
1789 : /* return x * n/d. x: rational; d,n,result: integral; d,n coprime */
1790 : static GEN
1791 158053 : Q_divmuli_to_int(GEN x, GEN d, GEN n)
1792 : {
1793 : long i, l;
1794 : GEN y, xn, xd;
1795 : pari_sp av;
1796 :
1797 158053 : switch(typ(x))
1798 : {
1799 : case t_INT:
1800 53746 : av = avma; y = diviiexact(x,d);
1801 53746 : return gerepileuptoint(av, mulii(y,n));
1802 :
1803 : case t_FRAC:
1804 90606 : xn = gel(x,1);
1805 90606 : xd = gel(x,2); av = avma;
1806 90606 : y = mulii(diviiexact(xn, d), diviiexact(n, xd));
1807 90606 : return gerepileuptoint(av, y);
1808 :
1809 : case t_VEC: case t_COL: case t_MAT:
1810 12035 : y = cgetg_copy(x, &l);
1811 12035 : for (i=1; i<l; i++) gel(y,i) = Q_divmuli_to_int(gel(x,i), d,n);
1812 12035 : return y;
1813 :
1814 : case t_POL:
1815 1645 : y = cgetg_copy(x, &l); y[1] = x[1];
1816 1645 : for (i=2; i<l; i++) gel(y,i) = Q_divmuli_to_int(gel(x,i), d,n);
1817 1645 : return y;
1818 :
1819 : case t_POLMOD:
1820 21 : y = cgetg(3, t_POLMOD);
1821 21 : gel(y,1) = RgX_copy(gel(x,1));
1822 21 : gel(y,2) = Q_divmuli_to_int(gel(x,2), d,n);
1823 21 : return y;
1824 : }
1825 0 : pari_err_TYPE("Q_divmuli_to_int",x);
1826 : return NULL; /* LCOV_EXCL_LINE */
1827 : }
1828 :
1829 : /* return x / d. x: rational; d,result: integral. */
1830 : static GEN
1831 11690672 : Q_divi_to_int(GEN x, GEN d)
1832 : {
1833 : long i, l;
1834 : GEN y;
1835 :
1836 11690672 : switch(typ(x))
1837 : {
1838 : case t_INT:
1839 9558882 : return diviiexact(x,d);
1840 :
1841 : case t_VEC: case t_COL: case t_MAT:
1842 744271 : y = cgetg_copy(x, &l);
1843 744271 : for (i=1; i<l; i++) gel(y,i) = Q_divi_to_int(gel(x,i), d);
1844 744271 : return y;
1845 :
1846 : case t_POL:
1847 1383312 : y = cgetg_copy(x, &l); y[1] = x[1];
1848 1383312 : for (i=2; i<l; i++) gel(y,i) = Q_divi_to_int(gel(x,i), d);
1849 1383312 : return y;
1850 :
1851 : case t_POLMOD:
1852 4207 : y = cgetg(3, t_POLMOD);
1853 4207 : gel(y,1) = RgX_copy(gel(x,1));
1854 4207 : gel(y,2) = Q_divi_to_int(gel(x,2), d);
1855 4207 : return y;
1856 : }
1857 0 : pari_err_TYPE("Q_divi_to_int",x);
1858 : return NULL; /* LCOV_EXCL_LINE */
1859 : }
1860 : /* c t_FRAC */
1861 : static GEN
1862 218702 : Q_divq_to_int(GEN x, GEN c)
1863 : {
1864 218702 : GEN n = gel(c,1), d = gel(c,2);
1865 218702 : if (is_pm1(n)) {
1866 211855 : GEN y = Q_muli_to_int(x,d);
1867 211855 : if (signe(n) < 0) y = gneg(y);
1868 211855 : return y;
1869 : }
1870 6847 : return Q_divmuli_to_int(x, n,d);
1871 : }
1872 :
1873 : /* return y = x / c, assuming x,c rational, and y integral */
1874 : GEN
1875 3073 : Q_div_to_int(GEN x, GEN c)
1876 : {
1877 3073 : switch(typ(c))
1878 : {
1879 1953 : case t_INT: return Q_divi_to_int(x, c);
1880 1120 : case t_FRAC: return Q_divq_to_int(x, c);
1881 : }
1882 0 : pari_err_TYPE("Q_div_to_int",c);
1883 : return NULL; /* LCOV_EXCL_LINE */
1884 : }
1885 : /* return y = x * c, assuming x,c rational, and y integral */
1886 : GEN
1887 0 : Q_mul_to_int(GEN x, GEN c)
1888 : {
1889 : GEN d, n;
1890 0 : switch(typ(c))
1891 : {
1892 0 : case t_INT: return Q_muli_to_int(x, c);
1893 : case t_FRAC:
1894 0 : n = gel(c,1);
1895 0 : d = gel(c,2);
1896 0 : return Q_divmuli_to_int(x, d,n);
1897 : }
1898 0 : pari_err_TYPE("Q_mul_to_int",c);
1899 : return NULL; /* LCOV_EXCL_LINE */
1900 : }
1901 :
1902 : GEN
1903 10300354 : Q_primitive_part(GEN x, GEN *ptc)
1904 : {
1905 10300354 : pari_sp av = avma;
1906 10300354 : GEN c = Q_content(x);
1907 10300354 : if (typ(c) == t_INT)
1908 : {
1909 10082772 : if (is_pm1(c)) { avma = av; c = NULL; }
1910 1596063 : else if (signe(c)) x = Q_divi_to_int(x, c);
1911 : }
1912 217582 : else x = Q_divq_to_int(x, c);
1913 10300354 : if (ptc) *ptc = c;
1914 10300354 : return x;
1915 : }
1916 : GEN
1917 1379778 : Q_primpart(GEN x) { return Q_primitive_part(x, NULL); }
1918 :
1919 : GEN
1920 98967 : vec_Q_primpart(GEN M)
1921 : {
1922 : long i, l;
1923 98967 : GEN N = cgetg_copy(M, &l);
1924 98967 : for (i = 1; i < l; i++) gel(N,i) = Q_primpart(gel(M,i));
1925 98967 : return N;
1926 : }
1927 :
1928 : /*******************************************************************/
1929 : /* */
1930 : /* SUBRESULTANT */
1931 : /* */
1932 : /*******************************************************************/
1933 : /* for internal use */
1934 : GEN
1935 2006634 : gdivexact(GEN x, GEN y)
1936 : {
1937 : long i,lx;
1938 : GEN z;
1939 2006634 : if (gequal1(y)) return x;
1940 1988099 : switch(typ(x))
1941 : {
1942 : case t_INT:
1943 1613966 : if (typ(y)==t_INT) return diviiexact(x,y);
1944 70 : if (!signe(x)) return gen_0;
1945 0 : break;
1946 : case t_INTMOD:
1947 14760 : case t_POLMOD: return gmul(x,ginv(y));
1948 : case t_POL:
1949 355397 : switch(typ(y))
1950 : {
1951 : case t_INTMOD:
1952 0 : case t_POLMOD: return gmul(x,ginv(y));
1953 : case t_POL: { /* not stack-clean */
1954 : long v;
1955 20706 : if (varn(x)!=varn(y)) break;
1956 19859 : v = RgX_valrem(y,&y);
1957 19859 : if (v) x = RgX_shift_shallow(x,-v);
1958 19859 : if (!degpol(y)) { y = gel(y,2); break; }
1959 16310 : return RgX_div(x,y);
1960 : }
1961 : }
1962 339087 : return RgX_Rg_divexact(x, y);
1963 : case t_VEC: case t_COL: case t_MAT:
1964 3976 : lx = lg(x); z = new_chunk(lx);
1965 3976 : for (i=1; i<lx; i++) gel(z,i) = gdivexact(gel(x,i),y);
1966 3976 : z[0] = x[0]; return z;
1967 : }
1968 0 : if (DEBUGLEVEL) pari_warn(warner,"missing case in gdivexact");
1969 0 : return gdiv(x,y);
1970 : }
1971 :
1972 : static GEN
1973 62863 : init_resultant(GEN x, GEN y)
1974 : {
1975 62863 : long tx = typ(x), ty = typ(y), vx, vy;
1976 62863 : if (is_scalar_t(tx) || is_scalar_t(ty))
1977 : {
1978 196 : if (gequal0(x) || gequal0(y)) return gmul(x,y); /* keep type info */
1979 196 : if (tx==t_POL) return gpowgs(y, degpol(x));
1980 0 : if (ty==t_POL) return gpowgs(x, degpol(y));
1981 0 : return gen_1;
1982 : }
1983 62667 : if (tx!=t_POL) pari_err_TYPE("resultant_all",x);
1984 62667 : if (ty!=t_POL) pari_err_TYPE("resultant_all",y);
1985 62667 : if (!signe(x) || !signe(y)) return gmul(RgX_get_0(x),RgX_get_0(y)); /*type*/
1986 62590 : vx = varn(x);
1987 62590 : vy = varn(y); if (vx == vy) return NULL;
1988 0 : return (varncmp(vx,vy) < 0)? gpowgs(y,degpol(x)): gpowgs(x,degpol(y));
1989 : }
1990 :
1991 : static long
1992 102969 : RgX_simpletype(GEN x)
1993 : {
1994 102969 : long T = t_INT, i, lx = lg(x);
1995 601343 : for (i = 2; i < lx; i++)
1996 : {
1997 507033 : GEN c = gel(x,i);
1998 507033 : long tc = typ(c);
1999 507033 : switch(tc) {
2000 : case t_INT:
2001 433747 : break;
2002 : case t_FRAC:
2003 33002 : if (T == t_INT) T = t_FRAC;
2004 33002 : break;
2005 : default:
2006 40284 : if (isinexact(c)) return t_REAL;
2007 31625 : T = 0; break;
2008 : }
2009 : }
2010 94310 : return T;
2011 : }
2012 :
2013 : /* x an RgX, y a scalar */
2014 : static GEN
2015 0 : scalar_res(GEN x, GEN y, GEN *U, GEN *V)
2016 : {
2017 0 : *V = gpowgs(y,degpol(x)-1);
2018 0 : *U = gen_0; return gmul(y, *V);
2019 : }
2020 :
2021 : /* return 0 if the subresultant chain can be interrupted.
2022 : * Set u = NULL if the resultant is 0. */
2023 : static int
2024 53698 : subres_step(GEN *u, GEN *v, GEN *g, GEN *h, GEN *uze, GEN *um1, long *signh)
2025 : {
2026 53698 : GEN u0, c, r, q = RgX_pseudodivrem(*u,*v, &r);
2027 : long du, dv, dr, degq;
2028 :
2029 53698 : if (gequal0(leading_coeff(r))) r = RgX_renormalize(r);
2030 53698 : dr = lg(r); if (!signe(r)) { *u = NULL; return 0; }
2031 53572 : du = degpol(*u);
2032 53572 : dv = degpol(*v);
2033 53572 : degq = du - dv;
2034 53572 : if (*um1 == gen_1)
2035 24935 : u0 = gpowgs(gel(*v,dv+2),degq+1);
2036 28637 : else if (*um1 == gen_0)
2037 13054 : u0 = gen_0;
2038 : else /* except in those 2 cases, um1 is an RgX */
2039 15583 : u0 = RgX_Rg_mul(*um1, gpowgs(gel(*v,dv+2),degq+1));
2040 :
2041 53572 : if (*uze == gen_0) /* except in that case, uze is an RgX */
2042 24935 : u0 = scalarpol(u0, varn(*u)); /* now an RgX */
2043 : else
2044 28637 : u0 = gsub(u0, gmul(q,*uze));
2045 :
2046 53572 : *um1 = *uze;
2047 53572 : *uze = u0; /* uze <- lead(v)^(degq + 1) * um1 - q * uze */
2048 :
2049 53572 : *u = *v; c = *g; *g = leading_coeff(*u);
2050 53572 : switch(degq)
2051 : {
2052 483 : case 0: break;
2053 : case 1:
2054 35716 : c = gmul(*h,c); *h = *g; break;
2055 : default:
2056 17373 : c = gmul(gpowgs(*h,degq), c);
2057 17373 : *h = gdivexact(gpowgs(*g,degq), gpowgs(*h,degq-1));
2058 : }
2059 53572 : *v = RgX_Rg_divexact(r,c);
2060 53572 : *uze= RgX_Rg_divexact(*uze,c);
2061 53572 : if (both_odd(du, dv)) *signh = -*signh;
2062 53572 : return (dr > 3);
2063 : }
2064 :
2065 : /* compute U, V s.t Ux + Vy = resultant(x,y) */
2066 : static GEN
2067 25644 : subresext_i(GEN x, GEN y, GEN *U, GEN *V)
2068 : {
2069 : pari_sp av, av2;
2070 25644 : long dx, dy, du, signh, tx = typ(x), ty = typ(y);
2071 : GEN r, z, g, h, p1, cu, cv, u, v, um1, uze, vze;
2072 :
2073 25644 : if (!is_extscalar_t(tx)) pari_err_TYPE("subresext",x);
2074 25644 : if (!is_extscalar_t(ty)) pari_err_TYPE("subresext",y);
2075 25644 : if (gequal0(x) || gequal0(y)) { *U = *V = gen_0; return gen_0; }
2076 25637 : if (tx != t_POL) {
2077 0 : if (ty != t_POL) { *U = ginv(x); *V = gen_0; return gen_1; }
2078 0 : return scalar_res(y,x,V,U);
2079 : }
2080 25637 : if (ty != t_POL) return scalar_res(x,y,U,V);
2081 25637 : if (varn(x) != varn(y))
2082 0 : return varncmp(varn(x), varn(y)) < 0? scalar_res(x,y,U,V)
2083 0 : : scalar_res(y,x,V,U);
2084 25637 : if (gequal0(leading_coeff(x))) x = RgX_renormalize(x);
2085 25637 : if (gequal0(leading_coeff(y))) y = RgX_renormalize(y);
2086 25637 : dx = degpol(x);
2087 25637 : dy = degpol(y);
2088 25637 : signh = 1;
2089 25637 : if (dx < dy)
2090 : {
2091 24223 : pswap(U,V); lswap(dx,dy); swap(x,y);
2092 24223 : if (both_odd(dx, dy)) signh = -signh;
2093 : }
2094 25637 : if (dy == 0)
2095 : {
2096 793 : *V = gpowgs(gel(y,2),dx-1);
2097 793 : *U = gen_0; return gmul(*V,gel(y,2));
2098 : }
2099 24844 : av = avma;
2100 24844 : u = x = primitive_part(x, &cu);
2101 24844 : v = y = primitive_part(y, &cv);
2102 24844 : g = h = gen_1; av2 = avma;
2103 24844 : um1 = gen_1; uze = gen_0;
2104 : for(;;)
2105 : {
2106 53404 : if (!subres_step(&u, &v, &g, &h, &uze, &um1, &signh)) break;
2107 28560 : if (gc_needed(av2,1))
2108 : {
2109 0 : if(DEBUGMEM>1) pari_warn(warnmem,"subresext, dr = %ld", degpol(v));
2110 0 : gerepileall(av2,6, &u,&v,&g,&h,&uze,&um1);
2111 : }
2112 28560 : }
2113 : /* uze an RgX */
2114 24844 : if (!u) { *U = *V = gen_0; avma = av; return gen_0; }
2115 24837 : z = gel(v,2); du = degpol(u);
2116 24837 : if (du > 1)
2117 : { /* z = gdivexact(gpowgs(z,du), gpowgs(h,du-1)); */
2118 7535 : p1 = gpowgs(gdiv(z,h),du-1);
2119 7535 : z = gmul(z,p1);
2120 7535 : uze = RgX_Rg_mul(uze, p1);
2121 : }
2122 24837 : if (signh < 0) { z = gneg_i(z); uze = RgX_neg(uze); }
2123 :
2124 24837 : vze = RgX_divrem(Rg_RgX_sub(z, RgX_mul(uze,x)), y, &r);
2125 24837 : if (signe(r)) pari_warn(warner,"inexact computation in subresext");
2126 : /* uze ppart(x) + vze ppart(y) = z = resultant(ppart(x), ppart(y)), */
2127 24837 : p1 = gen_1;
2128 24837 : if (cu) p1 = gmul(p1, gpowgs(cu,dy));
2129 24837 : if (cv) p1 = gmul(p1, gpowgs(cv,dx));
2130 24837 : cu = cu? gdiv(p1,cu): p1;
2131 24837 : cv = cv? gdiv(p1,cv): p1;
2132 24837 : z = gmul(z,p1);
2133 24837 : *U = RgX_Rg_mul(uze,cu);
2134 24837 : *V = RgX_Rg_mul(vze,cv);
2135 24837 : return z;
2136 : }
2137 : GEN
2138 0 : subresext(GEN x, GEN y, GEN *U, GEN *V)
2139 : {
2140 0 : pari_sp av = avma;
2141 0 : GEN z = subresext_i(x, y, U, V);
2142 0 : gerepileall(av, 3, &z, U, V);
2143 0 : return z;
2144 : }
2145 :
2146 : static GEN
2147 28 : zero_extgcd(GEN y, GEN *U, GEN *V, long vx)
2148 : {
2149 28 : GEN x=content(y);
2150 28 : *U=pol_0(vx); *V = scalarpol(ginv(x), vx); return gmul(y,*V);
2151 : }
2152 :
2153 : static int
2154 3437 : must_negate(GEN x)
2155 : {
2156 3437 : GEN t = leading_coeff(x);
2157 3437 : switch(typ(t))
2158 : {
2159 : case t_INT: case t_REAL:
2160 1141 : return (signe(t) < 0);
2161 : case t_FRAC:
2162 112 : return (signe(gel(t,1)) < 0);
2163 : }
2164 2184 : return 0;
2165 : }
2166 :
2167 : /* compute U, V s.t Ux + Vy = GCD(x,y) using subresultant */
2168 : GEN
2169 343 : RgX_extgcd(GEN x, GEN y, GEN *U, GEN *V)
2170 : {
2171 : pari_sp av, av2, tetpil;
2172 : long signh; /* junk */
2173 343 : long dx, dy, vx, tx = typ(x), ty = typ(y);
2174 : GEN z, g, h, p1, cu, cv, u, v, um1, uze, vze, *gptr[3];
2175 :
2176 343 : if (tx!=t_POL) pari_err_TYPE("RgX_extgcd",x);
2177 343 : if (ty!=t_POL) pari_err_TYPE("RgX_extgcd",y);
2178 343 : if ( varncmp(varn(x),varn(y))) pari_err_VAR("RgX_extgcd",x,y);
2179 343 : vx=varn(x);
2180 343 : if (!signe(x))
2181 : {
2182 14 : if (signe(y)) return zero_extgcd(y,U,V,vx);
2183 7 : *U = pol_0(vx); *V = pol_0(vx);
2184 7 : return pol_0(vx);
2185 : }
2186 329 : if (!signe(y)) return zero_extgcd(x,V,U,vx);
2187 308 : dx = degpol(x); dy = degpol(y);
2188 308 : if (dx < dy) { pswap(U,V); lswap(dx,dy); swap(x,y); }
2189 308 : if (dy==0) { *U=pol_0(vx); *V=ginv(y); return pol_1(vx); }
2190 :
2191 161 : av = avma;
2192 161 : u = x = primitive_part(x, &cu);
2193 161 : v = y = primitive_part(y, &cv);
2194 161 : g = h = gen_1; av2 = avma;
2195 161 : um1 = gen_1; uze = gen_0;
2196 : for(;;)
2197 : {
2198 182 : if (!subres_step(&u, &v, &g, &h, &uze, &um1, &signh)) break;
2199 21 : if (gc_needed(av2,1))
2200 : {
2201 0 : if(DEBUGMEM>1) pari_warn(warnmem,"RgX_extgcd, dr = %ld",degpol(v));
2202 0 : gerepileall(av2,6,&u,&v,&g,&h,&uze,&um1);
2203 : }
2204 21 : }
2205 161 : if (uze != gen_0) {
2206 : GEN r;
2207 42 : vze = RgX_divrem(RgX_sub(v, RgX_mul(uze,x)), y, &r);
2208 42 : if (signe(r)) pari_warn(warner,"inexact computation in RgX_extgcd");
2209 42 : if (cu) uze = RgX_Rg_div(uze,cu);
2210 42 : if (cv) vze = RgX_Rg_div(vze,cv);
2211 42 : p1 = ginv(content(v));
2212 : }
2213 : else /* y | x */
2214 : {
2215 119 : vze = cv ? RgX_Rg_div(pol_1(vx),cv): pol_1(vx);
2216 119 : uze = pol_0(vx);
2217 119 : p1 = gen_1;
2218 : }
2219 161 : if (must_negate(v)) p1 = gneg(p1);
2220 161 : tetpil = avma;
2221 161 : z = RgX_Rg_mul(v,p1);
2222 161 : *U = RgX_Rg_mul(uze,p1);
2223 161 : *V = RgX_Rg_mul(vze,p1);
2224 161 : gptr[0] = &z;
2225 161 : gptr[1] = U;
2226 161 : gptr[2] = V;
2227 161 : gerepilemanysp(av,tetpil,gptr,3); return z;
2228 : }
2229 :
2230 : int
2231 56 : RgXQ_ratlift(GEN x, GEN T, long amax, long bmax, GEN *P, GEN *Q)
2232 : {
2233 56 : pari_sp av = avma, av2, tetpil;
2234 : long signh; /* junk */
2235 : long vx;
2236 : GEN g, h, p1, cu, cv, u, v, um1, uze, *gptr[2];
2237 :
2238 56 : if (typ(x)!=t_POL) pari_err_TYPE("RgXQ_ratlift",x);
2239 56 : if (typ(T)!=t_POL) pari_err_TYPE("RgXQ_ratlift",T);
2240 56 : if ( varncmp(varn(x),varn(T)) ) pari_err_VAR("RgXQ_ratlift",x,T);
2241 56 : if (bmax < 0) pari_err_DOMAIN("ratlift", "bmax", "<", gen_0, stoi(bmax));
2242 56 : if (!signe(T)) {
2243 0 : if (degpol(x) <= amax) {
2244 0 : *P = RgX_copy(x);
2245 0 : *Q = pol_1(varn(x));
2246 0 : return 1;
2247 : }
2248 0 : return 0;
2249 : }
2250 56 : if (amax+bmax >= degpol(T))
2251 0 : pari_err_DOMAIN("ratlift", "amax+bmax", ">=", stoi(degpol(T)),
2252 : mkvec3(stoi(amax), stoi(bmax), T));
2253 56 : vx = varn(T);
2254 56 : u = x = primitive_part(x, &cu);
2255 56 : v = T = primitive_part(T, &cv);
2256 56 : g = h = gen_1; av2 = avma;
2257 56 : um1 = gen_1; uze = gen_0;
2258 : for(;;)
2259 : {
2260 112 : (void) subres_step(&u, &v, &g, &h, &uze, &um1, &signh);
2261 112 : if (!u || (typ(uze)==t_POL && degpol(uze)>bmax)) { avma=av; return 0; }
2262 112 : if (typ(v)!=t_POL || degpol(v)<=amax) break;
2263 56 : if (gc_needed(av2,1))
2264 : {
2265 0 : if(DEBUGMEM>1) pari_warn(warnmem,"RgXQ_ratlift, dr = %ld", degpol(v));
2266 0 : gerepileall(av2,6,&u,&v,&g,&h,&uze,&um1);
2267 : }
2268 56 : }
2269 56 : if (uze == gen_0)
2270 : {
2271 0 : avma = av; *P = pol_0(vx); *Q = pol_1(vx);
2272 0 : return 1;
2273 : }
2274 56 : if (cu) uze = RgX_Rg_div(uze,cu);
2275 56 : p1 = ginv(content(v));
2276 56 : if (must_negate(v)) p1 = gneg(p1);
2277 56 : tetpil = avma;
2278 56 : *P = RgX_Rg_mul(v,p1);
2279 56 : *Q = RgX_Rg_mul(uze,p1);
2280 56 : gptr[0] = P;
2281 56 : gptr[1] = Q;
2282 56 : gerepilemanysp(av,tetpil,gptr,2); return 1;
2283 : }
2284 :
2285 : /*******************************************************************/
2286 : /* */
2287 : /* RESULTANT USING DUCOS VARIANT */
2288 : /* */
2289 : /*******************************************************************/
2290 : /* x^n / y^(n-1), assume n > 0 */
2291 : static GEN
2292 20551 : Lazard(GEN x, GEN y, long n)
2293 : {
2294 : long a;
2295 : GEN c;
2296 :
2297 20551 : if (n == 1) return x;
2298 1327 : a = 1 << expu(n); /* a = 2^k <= n < 2^(k+1) */
2299 1327 : c=x; n-=a;
2300 4604 : while (a>1)
2301 : {
2302 1950 : a>>=1; c=gdivexact(gsqr(c),y);
2303 1950 : if (n>=a) { c=gdivexact(gmul(c,x),y); n -= a; }
2304 : }
2305 1327 : return c;
2306 : }
2307 :
2308 : /* F (x/y)^(n-1), assume n >= 1 */
2309 : static GEN
2310 20335 : Lazard2(GEN F, GEN x, GEN y, long n)
2311 : {
2312 20335 : if (n == 1) return F;
2313 707 : return RgX_Rg_divexact(RgX_Rg_mul(F, Lazard(x,y,n-1)), y);
2314 : }
2315 :
2316 : static GEN
2317 20335 : RgX_neg_i(GEN x, long lx)
2318 : {
2319 : long i;
2320 20335 : GEN y = cgetg(lx, t_POL); y[1] = x[1];
2321 20335 : for (i=2; i<lx; i++) gel(y,i) = gneg(gel(x,i));
2322 20335 : return y;
2323 : }
2324 : static GEN
2325 61306 : RgX_Rg_mul_i(GEN y, GEN x, long ly)
2326 : {
2327 : long i;
2328 : GEN z;
2329 61306 : if (isrationalzero(x)) return pol_0(varn(y));
2330 61292 : z = cgetg(ly,t_POL); z[1] = y[1];
2331 61292 : for (i = 2; i < ly; i++) gel(z,i) = gmul(x,gel(y,i));
2332 61292 : return z;
2333 : }
2334 : static long
2335 58429 : reductum_lg(GEN x, long lx)
2336 : {
2337 58429 : long i = lx-2;
2338 58429 : while (i > 1 && gequal0(gel(x,i))) i--;
2339 58429 : return i+1;
2340 : }
2341 :
2342 : #define addshift(x,y) RgX_addmulXn_shallow((x),(y),1)
2343 : /* delta = deg(P) - deg(Q) > 0, deg(Q) > 0, P,Q,Z t_POL in the same variable,
2344 : * s "scalar". Return prem(P, -Q) / s^delta lc(P) */
2345 : static GEN
2346 20335 : nextSousResultant(GEN P, GEN Q, GEN Z, GEN s)
2347 : {
2348 20335 : GEN p0, q0, h0, TMP, H, A, z0 = leading_coeff(Z);
2349 : long p, q, j, lP, lQ;
2350 : pari_sp av;
2351 :
2352 20335 : p = degpol(P); p0 = gel(P,p+2); lP = reductum_lg(P,lg(P));
2353 20335 : q = degpol(Q); q0 = gel(Q,q+2); lQ = reductum_lg(Q,lg(Q));
2354 : /* p > q. Very often p - 1 = q */
2355 20335 : av = avma;
2356 : /* H = RgX_neg(reductum(Z)) optimized, using Q ~ Z */
2357 20335 : H = RgX_neg_i(Z, lQ); /* deg H < q */
2358 :
2359 20335 : A = (q+2 < lP)? RgX_Rg_mul_i(H, gel(P,q+2), lQ): NULL;
2360 23947 : for (j = q+1; j < p; j++)
2361 : {
2362 3612 : if (degpol(H) == q-1)
2363 : { /* h0 = coeff of degree q-1 = leading coeff */
2364 3185 : h0 = gel(H,q+1); (void)normalizepol_lg(H, q+1);
2365 3185 : H = addshift(H, RgX_Rg_divexact(RgX_Rg_mul_i(Q, gneg(h0), lQ), q0));
2366 : }
2367 : else
2368 427 : H = RgX_shift_shallow(H, 1);
2369 3612 : if (j+2 < lP)
2370 : {
2371 3045 : TMP = RgX_Rg_mul(H, gel(P,j+2));
2372 3045 : A = A? RgX_add(A, TMP): TMP;
2373 : }
2374 3612 : if (gc_needed(av,1))
2375 : {
2376 147 : if(DEBUGMEM>1) pari_warn(warnmem,"nextSousResultant j = %ld/%ld",j,p);
2377 147 : gerepileall(av,A?2:1,&H,&A);
2378 : }
2379 : }
2380 20335 : if (q+2 < lP) lP = reductum_lg(P, q+3);
2381 20335 : TMP = RgX_Rg_mul_i(P, z0, lP);
2382 20335 : A = A? RgX_add(A, TMP): TMP;
2383 20335 : A = RgX_Rg_divexact(A, p0);
2384 20335 : if (degpol(H) == q-1)
2385 : {
2386 20027 : h0 = gel(H,q+1); (void)normalizepol_lg(H, q+1); /* destroy old H */
2387 20027 : A = RgX_add(RgX_Rg_mul(addshift(H,A),q0), RgX_Rg_mul_i(Q, gneg(h0), lQ));
2388 : }
2389 : else
2390 308 : A = RgX_Rg_mul(addshift(H,A), q0);
2391 20335 : return RgX_Rg_divexact(A, s);
2392 : }
2393 : #undef addshift
2394 :
2395 : /* Ducos's subresultant */
2396 : GEN
2397 21314 : RgX_resultant_all(GEN P, GEN Q, GEN *sol)
2398 : {
2399 : pari_sp av, av2;
2400 21314 : long dP, dQ, delta, sig = 1;
2401 : GEN cP, cQ, Z, s;
2402 :
2403 21314 : dP = degpol(P);
2404 21314 : dQ = degpol(Q); delta = dP - dQ;
2405 21314 : if (delta < 0)
2406 : {
2407 2499 : if (both_odd(dP, dQ)) sig = -1;
2408 2499 : swap(P,Q); lswap(dP, dQ); delta = -delta;
2409 : }
2410 21314 : if (sol) *sol = gen_0;
2411 21314 : av = avma;
2412 21314 : if (dQ <= 0)
2413 : {
2414 1470 : if (dQ < 0) return RgX_get_0(P);
2415 1470 : s = gpowgs(gel(Q,2), dP);
2416 1470 : if (sig == -1) s = gerepileupto(av, gneg(s));
2417 1470 : return s;
2418 : }
2419 19844 : P = primitive_part(P, &cP);
2420 19844 : Q = primitive_part(Q, &cQ);
2421 19844 : av2 = avma;
2422 19844 : s = gpowgs(leading_coeff(Q),delta);
2423 19844 : if (both_odd(dP, dQ)) sig = -sig;
2424 19844 : Z = Q;
2425 19844 : Q = RgX_pseudorem(P, Q);
2426 19844 : P = Z;
2427 60023 : while(degpol(Q) > 0)
2428 : {
2429 20335 : delta = degpol(P) - degpol(Q); /* > 0 */
2430 20335 : Z = Lazard2(Q, leading_coeff(Q), s, delta);
2431 20335 : if (both_odd(degpol(P), degpol(Q))) sig = -sig;
2432 20335 : Q = nextSousResultant(P, Q, Z, s);
2433 20335 : P = Z;
2434 20335 : if (gc_needed(av,1))
2435 : {
2436 13 : if(DEBUGMEM>1) pari_warn(warnmem,"resultant_all, degpol Q = %ld",degpol(Q));
2437 13 : gerepileall(av2,2,&P,&Q);
2438 : }
2439 20335 : s = leading_coeff(P);
2440 : }
2441 19844 : if (!signe(Q)) { avma = av; return RgX_get_0(Q); }
2442 19844 : s = Lazard(leading_coeff(Q), s, degpol(P));
2443 19844 : if (sig == -1) s = gneg(s);
2444 19844 : if (cP) s = gmul(s, gpowgs(cP,dQ));
2445 19844 : if (cQ) s = gmul(s, gpowgs(cQ,dP));
2446 19844 : if (sol) { *sol = P; gerepileall(av, 2, &s, sol); return s; }
2447 18276 : return gerepilecopy(av, s);
2448 : }
2449 : /* Return resultant(P,Q). If sol != NULL: set *sol to the last non-constant
2450 : * polynomial in the prs IF the sequence was computed, and gen_0 otherwise.
2451 : * Uses Sylvester's matrix if P or Q inexact, a modular algorithm if they
2452 : * are in Q[X], and Ducos/Lazard optimization of the subresultant algorithm
2453 : * in the "generic" case. */
2454 : GEN
2455 54197 : resultant_all(GEN P, GEN Q, GEN *sol)
2456 : {
2457 : long TP, TQ;
2458 : GEN s;
2459 :
2460 54197 : if (sol) *sol = gen_0;
2461 54197 : if ((s = init_resultant(P,Q))) return s;
2462 53924 : if ((TP = RgX_simpletype(P)) == t_REAL || (TQ = RgX_simpletype(Q)) == t_REAL)
2463 8645 : return resultant2(P,Q); /* inexact */
2464 45279 : if (TP && TQ) /* rational */
2465 : {
2466 24630 : if (TP == t_INT && TQ == t_INT) return ZX_resultant(P,Q);
2467 10530 : return QX_resultant(P,Q);
2468 : }
2469 20649 : return RgX_resultant_all(P, Q, sol);
2470 : }
2471 :
2472 : /*******************************************************************/
2473 : /* */
2474 : /* RESULTANT USING SYLVESTER MATRIX */
2475 : /* */
2476 : /*******************************************************************/
2477 : static GEN
2478 0 : _pol_0(void)
2479 : {
2480 0 : GEN x = cgetg(3,t_POL);
2481 0 : x[1] = 0;
2482 0 : gel(x,2) = gen_0; return x;
2483 : }
2484 :
2485 : static GEN
2486 72100 : sylvester_col(GEN x, long j, long d, long D)
2487 : {
2488 72100 : GEN c = cgetg(d+1,t_COL);
2489 : long i;
2490 72100 : for (i=1; i< j; i++) gel(c,i) = gen_0;
2491 72100 : for ( ; i<=D; i++) gel(c,i) = gel(x,D-i+2);
2492 72100 : for ( ; i<=d; i++) gel(c,i) = gen_0;
2493 72100 : return c;
2494 : }
2495 :
2496 : GEN
2497 8673 : sylvestermatrix_i(GEN x, GEN y)
2498 : {
2499 : long j,d,dx,dy;
2500 : GEN M;
2501 :
2502 8673 : dx = degpol(x); if (dx < 0) { dx = 0; x = _pol_0(); }
2503 8673 : dy = degpol(y); if (dy < 0) { dy = 0; y = _pol_0(); }
2504 8673 : d = dx+dy; M = cgetg(d+1,t_MAT);
2505 8673 : for (j=1; j<=dy; j++) gel(M,j) = sylvester_col(x,j,d,j+dx);
2506 8673 : for (j=1; j<=dx; j++) gel(M,j+dy) = sylvester_col(y,j,d,j+dy);
2507 8673 : return M;
2508 : }
2509 :
2510 : GEN
2511 7 : sylvestermatrix(GEN x, GEN y)
2512 : {
2513 : long i,j,lx;
2514 7 : if (typ(x)!=t_POL) pari_err_TYPE("sylvestermatrix",x);
2515 7 : if (typ(y)!=t_POL) pari_err_TYPE("sylvestermatrix",y);
2516 7 : if (varn(x) != varn(y)) pari_err_VAR("sylvestermatrix",x,y);
2517 7 : x = sylvestermatrix_i(x,y); lx = lg(x);
2518 28 : for (i=1; i<lx; i++)
2519 21 : for (j=1; j<lx; j++) gcoeff(x,i,j) = gcopy(gcoeff(x,i,j));
2520 7 : return x;
2521 : }
2522 :
2523 : GEN
2524 8666 : resultant2(GEN x, GEN y)
2525 : {
2526 : pari_sp av;
2527 : GEN r;
2528 8666 : if ((r = init_resultant(x,y))) return r;
2529 8666 : av = avma; return gerepileupto(av,det(sylvestermatrix_i(x,y)));
2530 : }
2531 :
2532 : /* If x a t_POL, let vx = main variable of x; return a t_POL in variable v0:
2533 : * if vx <= v, return subst(x, v, pol_x(v0))
2534 : * if vx > v, return scalarpol(x, v0) */
2535 : static GEN
2536 2394 : fix_pol(GEN x, long v, long v0)
2537 : {
2538 : long vx;
2539 2394 : if (typ(x) != t_POL) return x;
2540 2394 : vx = varn(x);
2541 2394 : if (v == vx)
2542 : {
2543 2331 : if (v0 != v) { x = leafcopy(x); setvarn(x, v0); }
2544 2331 : return x;
2545 : }
2546 63 : if (varncmp(v, vx) > 0)
2547 : {
2548 56 : x = gsubst(x,v,pol_x(v0));
2549 56 : if (typ(x) == t_POL && varn(x) == v0) return x;
2550 : }
2551 35 : return scalarpol_shallow(x, v0);
2552 : }
2553 :
2554 : /* resultant of x and y with respect to variable v, or with respect to their
2555 : * main variable if v < 0. */
2556 : GEN
2557 1393 : polresultant0(GEN x, GEN y, long v, long flag)
2558 : {
2559 1393 : long v0 = 0;
2560 1393 : pari_sp av = avma;
2561 :
2562 1393 : if (v >= 0)
2563 : {
2564 1183 : v0 = fetch_var_higher();
2565 1183 : x = fix_pol(x,v, v0);
2566 1183 : y = fix_pol(y,v, v0);
2567 : }
2568 1393 : switch(flag)
2569 : {
2570 : case 2:
2571 1386 : case 0: x=resultant_all(x,y,NULL); break;
2572 7 : case 1: x=resultant2(x,y); break;
2573 0 : default: pari_err_FLAG("polresultant");
2574 : }
2575 1393 : if (v >= 0) (void)delete_var();
2576 1393 : return gerepileupto(av,x);
2577 : }
2578 : GEN
2579 868 : polresultantext0(GEN x, GEN y, long v)
2580 : {
2581 : GEN R, U, V;
2582 868 : long v0 = 0;
2583 868 : pari_sp av = avma;
2584 :
2585 868 : if (v >= 0)
2586 : {
2587 14 : v0 = fetch_var_higher();
2588 14 : x = fix_pol(x,v, v0);
2589 14 : y = fix_pol(y,v, v0);
2590 : }
2591 868 : R = subresext_i(x,y, &U,&V);
2592 868 : if (v >= 0)
2593 : {
2594 14 : (void)delete_var();
2595 14 : if (typ(U) == t_POL && varn(U) != v) U = poleval(U, pol_x(v));
2596 14 : if (typ(V) == t_POL && varn(V) != v) V = poleval(V, pol_x(v));
2597 : }
2598 868 : return gerepilecopy(av, mkvec3(U,V,R));
2599 : }
2600 : GEN
2601 840 : polresultantext(GEN x, GEN y) { return polresultantext0(x,y,-1); }
2602 :
2603 : /*******************************************************************/
2604 : /* */
2605 : /* CHARACTERISTIC POLYNOMIAL USING RESULTANT */
2606 : /* */
2607 : /*******************************************************************/
2608 :
2609 : /* (v - x)^d */
2610 : static GEN
2611 140 : caract_const(pari_sp av, GEN x, long v, long d)
2612 140 : { return gerepileupto(av, gpowgs(deg1pol_shallow(gen_1, gneg_i(x), v), d)); }
2613 :
2614 : /* return caract(Mod(x,T)) in variable v */
2615 : GEN
2616 14391 : RgXQ_charpoly(GEN x, GEN T, long v)
2617 : {
2618 14391 : pari_sp av = avma;
2619 14391 : long d = degpol(T), dx, vx, vp, v0;
2620 : GEN ch, L;
2621 :
2622 14391 : if (typ(x) != t_POL) return caract_const(av, x, v, d);
2623 14377 : vx = varn(x);
2624 14377 : vp = varn(T);
2625 14377 : if (varncmp(vx, vp) > 0) return caract_const(av, x, v, d);
2626 14377 : if (varncmp(vx, vp) < 0) pari_err_PRIORITY("RgXQ_charpoly", x, "<", vp);
2627 14377 : dx = degpol(x);
2628 14377 : if (dx >= degpol(T)) { x = RgX_rem(x, T); dx = degpol(x); }
2629 14377 : if (dx <= 0) return dx? pol_xn(d, v): caract_const(av, gel(x,2), v, d);
2630 :
2631 14321 : v0 = fetch_var_higher();
2632 14321 : x = RgX_neg(x);
2633 14321 : gel(x,2) = gadd(gel(x,2), pol_x(v));
2634 14321 : setvarn(x, v0);
2635 14321 : T = leafcopy(T); setvarn(T, v0);
2636 14321 : ch = resultant_all(T, x, NULL);
2637 14321 : (void)delete_var();
2638 : /* test for silly input: x mod (deg 0 polynomial) */
2639 14321 : if (typ(ch) != t_POL)
2640 7 : pari_err_PRIORITY("RgXQ_charpoly", pol_x(v), "<", gvar(ch));
2641 14314 : L = leading_coeff(ch);
2642 14314 : if (!gequal1(L)) ch = RgX_Rg_div(ch, L);
2643 14314 : return gerepileupto(av, ch);
2644 : }
2645 :
2646 : /* characteristic polynomial (in v) of x over nf, where x is an element of the
2647 : * algebra nf[t]/(Q(t)) */
2648 : GEN
2649 224 : rnfcharpoly(GEN nf, GEN Q, GEN x, long v)
2650 : {
2651 224 : const char *f = "rnfcharpoly";
2652 224 : long dQ = degpol(Q);
2653 224 : pari_sp av = avma;
2654 : GEN T;
2655 :
2656 224 : if (v < 0) v = 0;
2657 224 : nf = checknf(nf); T = nf_get_pol(nf);
2658 224 : Q = RgX_nffix(f, T,Q,0);
2659 224 : switch(typ(x))
2660 : {
2661 : case t_INT:
2662 28 : case t_FRAC: return caract_const(av, x, v, dQ);
2663 : case t_POLMOD:
2664 91 : x = polmod_nffix2(f,T,Q, x,0);
2665 49 : break;
2666 : case t_POL:
2667 56 : x = varn(x) == varn(T)? Rg_nffix(f,T,x,0): RgX_nffix(f, T,x,0);
2668 42 : break;
2669 49 : default: pari_err_TYPE(f,x);
2670 : }
2671 91 : if (typ(x) != t_POL) return caract_const(av, x, v, dQ);
2672 : /* x a t_POL in variable vQ */
2673 49 : if (degpol(x) >= dQ) x = RgX_rem(x, Q);
2674 49 : if (dQ <= 1) return caract_const(av, constant_coeff(x), v, 1);
2675 49 : return gerepilecopy(av, lift_if_rational( RgXQ_charpoly(x, Q, v) ));
2676 : }
2677 :
2678 : /*******************************************************************/
2679 : /* */
2680 : /* GCD USING SUBRESULTANT */
2681 : /* */
2682 : /*******************************************************************/
2683 : static int inexact(GEN x, int *simple, int *rational);
2684 : static int
2685 6764544 : isinexactall(GEN x, int *simple, int *rational)
2686 : {
2687 6764544 : long i, lx = lg(x);
2688 31851318 : for (i=2; i<lx; i++)
2689 25086774 : if (inexact(gel(x,i), simple, rational)) return 1;
2690 6764544 : return 0;
2691 : }
2692 : /* return 1 if coeff explosion is not possible */
2693 : static int
2694 25086830 : inexact(GEN x, int *simple, int *rational)
2695 : {
2696 25086830 : int junk = 0;
2697 25086830 : switch(typ(x))
2698 : {
2699 25045663 : case t_INT: case t_FRAC: return 0;
2700 :
2701 0 : case t_REAL: case t_PADIC: case t_SER: return 1;
2702 :
2703 : case t_INTMOD:
2704 : case t_FFELT:
2705 7287 : *rational = 0;
2706 7287 : if (!*simple) *simple = 1;
2707 7287 : return 0;
2708 :
2709 : case t_COMPLEX:
2710 0 : *rational = 0;
2711 0 : return inexact(gel(x,1), simple, rational)
2712 0 : || inexact(gel(x,2), simple, rational);
2713 : case t_QUAD:
2714 0 : *rational = *simple = 0;
2715 0 : return inexact(gel(x,2), &junk, rational)
2716 0 : || inexact(gel(x,3), &junk, rational);
2717 :
2718 : case t_POLMOD:
2719 4655 : *rational = 0;
2720 4655 : return isinexactall(gel(x,1), simple, rational);
2721 : case t_POL:
2722 29197 : *rational = 0;
2723 29197 : *simple = -1;
2724 29197 : return isinexactall(x, &junk, rational);
2725 : case t_RFRAC:
2726 28 : *rational = 0;
2727 28 : *simple = -1;
2728 56 : return inexact(gel(x,1), &junk, rational)
2729 28 : || inexact(gel(x,2), &junk, rational);
2730 : }
2731 0 : *rational = 0;
2732 0 : *simple = -1; return 0;
2733 : }
2734 :
2735 : /* x monomial, y t_POL in the same variable */
2736 : static GEN
2737 7426802 : gcdmonome(GEN x, GEN y)
2738 : {
2739 7426802 : pari_sp av = avma;
2740 7426802 : long dx = degpol(x), e = RgX_valrem(y, &y);
2741 7426802 : long i, l = lg(y);
2742 7426802 : GEN t, v = cgetg(l, t_VEC);
2743 7426802 : gel(v,1) = gel(x,dx+2);
2744 7426802 : for (i = 2; i < l; i++) gel(v,i) = gel(y,i);
2745 7426802 : t = content(v); /* gcd(lc(x), cont(y)) */
2746 7426802 : t = simplify_shallow(t);
2747 7426802 : if (dx < e) e = dx;
2748 7426802 : return gerepileupto(av, monomialcopy(t, e, varn(x)));
2749 : }
2750 :
2751 : /* x, y are t_POL in the same variable */
2752 : GEN
2753 10792519 : RgX_gcd(GEN x, GEN y)
2754 : {
2755 : long dx, dy;
2756 : pari_sp av, av1;
2757 : GEN d, g, h, p1, p2, u, v;
2758 10792519 : int simple = 0, rational = 1;
2759 :
2760 10792519 : if (isexactzero(y)) return RgX_copy(x);
2761 10792155 : if (isexactzero(x)) return RgX_copy(y);
2762 10792148 : if (RgX_is_monomial(x)) return gcdmonome(x,y);
2763 3995271 : if (RgX_is_monomial(y)) return gcdmonome(y,x);
2764 3365346 : if (isinexactall(x,&simple,&rational) || isinexactall(y,&simple,&rational))
2765 : {
2766 0 : av = avma; u = ggcd(content(x), content(y));
2767 0 : return gerepileupto(av, scalarpol(u, varn(x)));
2768 : }
2769 3365346 : if (rational) return QX_gcd(x,y); /* Q[X] */
2770 :
2771 4396 : av = avma;
2772 4396 : if (simple > 0) x = RgX_gcd_simple(x,y);
2773 : else
2774 : {
2775 3619 : dx = lg(x); dy = lg(y);
2776 3619 : if (dx < dy) { swap(x,y); lswap(dx,dy); }
2777 3619 : if (dy==3)
2778 : {
2779 0 : d = ggcd(gel(y,2), content(x));
2780 0 : return gerepileupto(av, scalarpol(d, varn(x)));
2781 : }
2782 3619 : u = primitive_part(x, &p1); if (!p1) p1 = gen_1;
2783 3619 : v = primitive_part(y, &p2); if (!p2) p2 = gen_1;
2784 3619 : d = ggcd(p1,p2);
2785 3619 : av1 = avma;
2786 3619 : g = h = gen_1;
2787 : for(;;)
2788 : {
2789 4753 : GEN r = RgX_pseudorem(u,v);
2790 4753 : long degq, du, dv, dr = lg(r);
2791 :
2792 4753 : if (!signe(r)) break;
2793 2310 : if (dr <= 3)
2794 : {
2795 1176 : avma = av1; return gerepileupto(av, scalarpol(d, varn(x)));
2796 : }
2797 1134 : if (DEBUGLEVEL > 9) err_printf("RgX_gcd: dr = %ld\n", degpol(r));
2798 1134 : du = lg(u); dv = lg(v); degq = du-dv;
2799 1134 : u = v; p1 = g; g = leading_coeff(u);
2800 1134 : switch(degq)
2801 : {
2802 189 : case 0: break;
2803 : case 1:
2804 763 : p1 = gmul(h,p1); h = g; break;
2805 : default:
2806 182 : p1 = gmul(gpowgs(h,degq), p1);
2807 182 : h = gdiv(gpowgs(g,degq), gpowgs(h,degq-1));
2808 : }
2809 1134 : v = RgX_Rg_div(r,p1);
2810 1134 : if (gc_needed(av1,1))
2811 : {
2812 0 : if(DEBUGMEM>1) pari_warn(warnmem,"RgX_gcd");
2813 0 : gerepileall(av1,4, &u,&v,&g,&h);
2814 : }
2815 1134 : }
2816 2443 : x = RgX_Rg_mul(primpart(v), d);
2817 : }
2818 3220 : if (must_negate(x)) x = RgX_neg(x);
2819 3220 : return gerepileupto(av,x);
2820 : }
2821 :
2822 : /* disc P = (-1)^(n(n-1)/2) lc(P)^(n - deg P' - 2) Res(P,P'), n = deg P */
2823 : static GEN
2824 4746 : RgX_disc_aux(GEN P)
2825 : {
2826 4746 : long n = degpol(P), TP, dd;
2827 : GEN D, L, y, p;
2828 4746 : if (!signe(P) || !n) return RgX_get_0(P);
2829 4746 : if (n == 1) return RgX_get_1(P);
2830 4634 : if (n == 2) {
2831 868 : GEN a = gel(P,4), b = gel(P,3), c = gel(P,2);
2832 868 : return gsub(gsqr(b), gmul2n(gmul(a,c),2));
2833 : }
2834 3766 : TP = RgX_simpletype(P);
2835 3766 : if (TP == t_INT) return ZX_disc(P);
2836 469 : if (TP == t_FRAC) return QX_disc(P);
2837 469 : p = NULL;
2838 469 : if (RgX_is_FpX(P, &p) && p)
2839 28 : return Fp_to_mod(FpX_disc(RgX_to_FpX(P,p), p), p);
2840 :
2841 441 : y = RgX_deriv(P);
2842 441 : if (!signe(y)) return RgX_get_0(y);
2843 441 : dd = degpol(P)-2 - degpol(y);
2844 441 : if (TP == t_REAL)
2845 14 : D = resultant2(P,y);
2846 : else
2847 : {
2848 427 : D = RgX_resultant_all(P, y, NULL);
2849 427 : if (D == gen_0) return RgX_get_0(y);
2850 : }
2851 441 : L = leading_coeff(P);
2852 441 : if (dd && !gequal1(L)) D = (dd == -1)? gdiv(D, L): gmul(D, gpowgs(L, dd));
2853 441 : if (n & 2) D = gneg(D);
2854 441 : return D;
2855 : }
2856 : GEN
2857 1687 : RgX_disc(GEN x) { pari_sp av = avma; return gerepileupto(av, RgX_disc_aux(x)); }
2858 :
2859 : GEN
2860 3073 : poldisc0(GEN x, long v)
2861 : {
2862 : pari_sp av;
2863 3073 : switch(typ(x))
2864 : {
2865 : case t_POL:
2866 : {
2867 : GEN D;
2868 3059 : long v0 = -1;
2869 3059 : av = avma;
2870 3059 : if (v >= 0 && v != varn(x))
2871 : {
2872 0 : v0 = fetch_var_higher();
2873 0 : x = fix_pol(x,v, v0);
2874 : }
2875 3059 : D = RgX_disc_aux(x);
2876 3059 : if (v0 >= 0) (void)delete_var();
2877 3059 : return gerepileupto(av, D);
2878 : }
2879 :
2880 : case t_COMPLEX:
2881 0 : return utoineg(4);
2882 :
2883 : case t_QUAD:
2884 0 : return quad_disc(x);
2885 : case t_POLMOD:
2886 0 : return poldisc0(gel(x,1), v);
2887 :
2888 : case t_QFR: case t_QFI:
2889 14 : av = avma; return gerepileuptoint(av, qfb_disc(x));
2890 :
2891 : case t_VEC: case t_COL: case t_MAT:
2892 : {
2893 : long i;
2894 0 : GEN z = cgetg_copy(x, &i);
2895 0 : for (i--; i; i--) gel(z,i) = poldisc0(gel(x,i), v);
2896 0 : return z;
2897 : }
2898 : }
2899 0 : pari_err_TYPE("poldisc",x);
2900 : return NULL; /* LCOV_EXCL_LINE */
2901 : }
2902 :
2903 : GEN
2904 7 : reduceddiscsmith(GEN x)
2905 : {
2906 7 : long j, n = degpol(x);
2907 7 : pari_sp av = avma;
2908 : GEN xp, M;
2909 :
2910 7 : if (typ(x) != t_POL) pari_err_TYPE("poldiscreduced",x);
2911 7 : if (n<=0) pari_err_CONSTPOL("poldiscreduced");
2912 7 : RgX_check_ZX(x,"poldiscreduced");
2913 7 : if (!gequal1(gel(x,n+2)))
2914 0 : pari_err_IMPL("non-monic polynomial in poldiscreduced");
2915 7 : M = cgetg(n+1,t_MAT);
2916 7 : xp = ZX_deriv(x);
2917 28 : for (j=1; j<=n; j++)
2918 : {
2919 21 : gel(M,j) = RgX_to_RgC(xp, n);
2920 21 : if (j<n) xp = RgX_rem(RgX_shift_shallow(xp, 1), x);
2921 : }
2922 7 : return gerepileupto(av, ZM_snf(M));
2923 : }
2924 :
2925 : /***********************************************************************/
2926 : /** **/
2927 : /** STURM ALGORITHM **/
2928 : /** (number of real roots of x in [a,b]) **/
2929 : /** **/
2930 : /***********************************************************************/
2931 : static GEN
2932 819 : R_to_Q_up(GEN x)
2933 : {
2934 : long e;
2935 819 : switch(typ(x))
2936 : {
2937 819 : case t_INT: case t_FRAC: case t_INFINITY: return x;
2938 : case t_REAL:
2939 0 : x = mantissa_real(x,&e);
2940 0 : return gmul2n(addiu(x,1), -e);
2941 0 : default: pari_err_TYPE("R_to_Q_up", x);
2942 : return NULL; /* LCOV_EXCL_LINE */
2943 : }
2944 : }
2945 : static GEN
2946 819 : R_to_Q_down(GEN x)
2947 : {
2948 : long e;
2949 819 : switch(typ(x))
2950 : {
2951 805 : case t_INT: case t_FRAC: case t_INFINITY: return x;
2952 : case t_REAL:
2953 14 : x = mantissa_real(x,&e);
2954 14 : return gmul2n(subiu(x,1), -e);
2955 0 : default: pari_err_TYPE("R_to_Q_down", x);
2956 : return NULL; /* LCOV_EXCL_LINE */
2957 : }
2958 : }
2959 :
2960 : static long
2961 833 : sturmpart_i(GEN x, GEN ab)
2962 : {
2963 833 : long tx = typ(x);
2964 833 : if (gequal0(x)) pari_err_ROOTS0("sturm");
2965 833 : if (tx != t_POL)
2966 : {
2967 0 : if (is_real_t(tx)) return 0;
2968 0 : pari_err_TYPE("sturm",x);
2969 : }
2970 833 : if (lg(x) == 3) return 0;
2971 833 : if (!RgX_is_ZX(x)) x = RgX_rescale_to_int(x);
2972 833 : if (!ZX_is_squarefree(x))
2973 14 : pari_err_DOMAIN("polsturm","issquarefree(pol)","=",gen_0,x);
2974 819 : if (ab)
2975 : {
2976 : GEN A, B;
2977 819 : if (typ(ab) != t_VEC || lg(ab) != 3) pari_err_TYPE("RgX_sturmpart", ab);
2978 819 : A = R_to_Q_down(gel(ab,1));
2979 819 : B = R_to_Q_up(gel(ab,2));
2980 819 : ab = mkvec2(A, B);
2981 : }
2982 819 : return ZX_sturmpart(x, ab);
2983 : }
2984 : /* Deprecated: RgX_sturmpart() should be preferred */
2985 : long
2986 833 : sturmpart(GEN x, GEN a, GEN b)
2987 : {
2988 833 : pari_sp av = avma;
2989 : long r;
2990 833 : if (!b && a && typ(a) == t_VEC) return RgX_sturmpart(x, a);
2991 700 : if (!a) a = mkmoo();
2992 700 : if (!b) b = mkoo();
2993 700 : r = sturmpart_i(x, mkvec2(a, b));
2994 686 : avma = av; return r;
2995 : }
2996 : long
2997 133 : RgX_sturmpart(GEN x, GEN ab)
2998 : {
2999 133 : pari_sp av = avma;
3000 133 : long r = sturmpart_i(x, ab);
3001 133 : avma = av; return r;
3002 : }
3003 :
3004 : /***********************************************************************/
3005 : /** **/
3006 : /** GENERIC EXTENDED GCD **/
3007 : /** **/
3008 : /***********************************************************************/
3009 : /* assume typ(x) = typ(y) = t_POL */
3010 : GEN
3011 24776 : RgXQ_inv(GEN x, GEN y)
3012 : {
3013 24776 : long vx=varn(x), vy=varn(y);
3014 : pari_sp av;
3015 : GEN u, v, d;
3016 :
3017 49552 : while (vx != vy)
3018 : {
3019 0 : if (varncmp(vx,vy) > 0)
3020 : {
3021 0 : d = (vx == NO_VARIABLE)? ginv(x): gred_rfrac_simple(gen_1, x);
3022 0 : return scalarpol(d, vy);
3023 : }
3024 0 : if (lg(x)!=3) pari_err_INV("RgXQ_inv",mkpolmod(x,y));
3025 0 : x = gel(x,2); vx = gvar(x);
3026 : }
3027 24776 : av = avma; d = subresext_i(x,y,&u,&v/*junk*/);
3028 24776 : if (gequal0(d)) pari_err_INV("RgXQ_inv",mkpolmod(x,y));
3029 24769 : d = gdiv(u,d);
3030 24769 : if (typ(d) != t_POL || varn(d) != vy) d = scalarpol(d, vy);
3031 24769 : return gerepileupto(av, d);
3032 : }
3033 :
3034 : /*Assume x is a polynomial and y is not */
3035 : static GEN
3036 112 : scalar_bezout(GEN x, GEN y, GEN *U, GEN *V)
3037 : {
3038 112 : long vx = varn(x);
3039 112 : int xis0 = signe(x)==0, yis0 = gequal0(y);
3040 112 : if (xis0 && yis0) { *U = *V = pol_0(vx); return pol_0(vx); }
3041 84 : if (yis0) { *U=pol_1(vx); *V = pol_0(vx); return RgX_copy(x);}
3042 56 : *U=pol_0(vx); *V= ginv(y); return pol_1(vx);
3043 : }
3044 : /* Assume x==0, y!=0 */
3045 : static GEN
3046 63 : zero_bezout(GEN y, GEN *U, GEN *V)
3047 : {
3048 63 : *U=gen_0; *V = ginv(y); return gen_1;
3049 : }
3050 :
3051 : GEN
3052 280 : gbezout(GEN x, GEN y, GEN *u, GEN *v)
3053 : {
3054 280 : long tx=typ(x), ty=typ(y), vx;
3055 280 : if (tx == t_INT && ty == t_INT) return bezout(x,y,u,v);
3056 245 : if (tx != t_POL)
3057 : {
3058 140 : if (ty == t_POL)
3059 56 : return scalar_bezout(y,x,v,u);
3060 : else
3061 : {
3062 84 : int xis0 = gequal0(x), yis0 = gequal0(y);
3063 84 : if (xis0 && yis0) { *u = *v = gen_0; return gen_0; }
3064 63 : if (xis0) return zero_bezout(y,u,v);
3065 42 : else return zero_bezout(x,v,u);
3066 : }
3067 : }
3068 105 : else if (ty != t_POL) return scalar_bezout(x,y,u,v);
3069 49 : vx = varn(x);
3070 49 : if (vx != varn(y))
3071 0 : return varncmp(vx, varn(y)) < 0? scalar_bezout(x,y,u,v)
3072 0 : : scalar_bezout(y,x,v,u);
3073 49 : return RgX_extgcd(x,y,u,v);
3074 : }
3075 :
3076 : GEN
3077 280 : gcdext0(GEN x, GEN y)
3078 : {
3079 280 : GEN z=cgetg(4,t_VEC);
3080 280 : gel(z,3) = gbezout(x,y,(GEN*)(z+1),(GEN*)(z+2));
3081 280 : return z;
3082 : }
3083 :
3084 : /*******************************************************************/
3085 : /* */
3086 : /* GENERIC (modular) INVERSE */
3087 : /* */
3088 : /*******************************************************************/
3089 :
3090 : GEN
3091 2660 : ginvmod(GEN x, GEN y)
3092 : {
3093 2660 : long tx=typ(x);
3094 :
3095 2660 : switch(typ(y))
3096 : {
3097 : case t_POL:
3098 2660 : if (tx==t_POL) return RgXQ_inv(x,y);
3099 1105 : if (is_scalar_t(tx)) return ginv(x);
3100 0 : break;
3101 : case t_INT:
3102 0 : if (tx==t_INT) return Fp_inv(x,y);
3103 0 : if (tx==t_POL) return gen_0;
3104 : }
3105 0 : pari_err_TYPE2("ginvmod",x,y);
3106 : return NULL; /* LCOV_EXCL_LINE */
3107 : }
3108 :
3109 : /***********************************************************************/
3110 : /** **/
3111 : /** NEWTON POLYGON **/
3112 : /** **/
3113 : /***********************************************************************/
3114 :
3115 : /* assume leading coeff of x is non-zero */
3116 : GEN
3117 28 : newtonpoly(GEN x, GEN p)
3118 : {
3119 : GEN y;
3120 : long n,ind,a,b,c,u1,u2,r1,r2;
3121 28 : long *vval, num[] = {evaltyp(t_INT) | _evallg(3), 0, 0};
3122 :
3123 28 : if (typ(x)!=t_POL) pari_err_TYPE("newtonpoly",x);
3124 28 : n=degpol(x); if (n<=0) return cgetg(1,t_VEC);
3125 28 : y = cgetg(n+1,t_VEC); x += 2; /* now x[i] = term of degree i */
3126 28 : vval = (long *) pari_malloc(sizeof(long)*(n+1));
3127 28 : for (a=0; a<=n; a++) vval[a] = gvaluation(gel(x,a),p);
3128 42 : for (a=0, ind=1; a<n; a++)
3129 : {
3130 42 : if (vval[a] != LONG_MAX) break;
3131 14 : gel(y,ind++) = mkoo();
3132 : }
3133 84 : for (b=a+1; b<=n; a=b, b=a+1)
3134 : {
3135 56 : while (vval[b] == LONG_MAX) b++;
3136 56 : u1 = vval[a]-vval[b];
3137 56 : u2 = b-a;
3138 154 : for (c=b+1; c<=n; c++)
3139 : {
3140 98 : if (vval[c] == LONG_MAX) continue;
3141 70 : r1 = vval[a]-vval[c];
3142 70 : r2 = c-a;
3143 70 : if (u1*r2 <= u2*r1) { u1 = r1; u2 = r2; b = c; }
3144 : }
3145 56 : while (ind<=b) { affsi(u1,num); gel(y,ind++) = gdivgs(num,u2); }
3146 : }
3147 28 : pari_free(vval); return y;
3148 : }
|
