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