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