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