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