1 /* $Id: polarit2.c 8460 2007-03-27 15:53:51Z kb $
2
3 Copyright (C) 2000 The PARI group.
4
5 This file is part of the PARI/GP package.
6
7 PARI/GP is free software; you can redistribute it and/or modify it under the
8 terms of the GNU General Public License as published by the Free Software
9 Foundation. It is distributed in the hope that it will be useful, but WITHOUT
10 ANY WARRANTY WHATSOEVER.
11
12 Check the License for details. You should have received a copy of it, along
13 with the package; see the file 'COPYING'. If not, write to the Free Software
14 Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. */
15
16 /***********************************************************************/
17 /** **/
18 /** ARITHMETIC OPERATIONS ON POLYNOMIALS **/
19 /** (second part) **/
20 /** **/
21 /***********************************************************************/
22 #include "pari.h"
23 #include "paripriv.h"
24
25 #define addshift(x,y) addshiftpol((x),(y),1)
26
27 /* compute Newton sums S_1(P), ... , S_n(P). S_k(P) = sum a_j^k, a_j root of P
28 * If N != NULL, assume p-adic roots and compute mod N [assume integer coeffs]
29 * If T != NULL, compute mod (T,N) [assume integer coeffs if N != NULL]
30 * If y0!= NULL, precomputed i-th powers, i=1..m, m = length(y0).
31 * Not memory clean in the latter case */
32 GEN
polsym_gen(GEN P,GEN y0,long n,GEN T,GEN N)33 polsym_gen(GEN P, GEN y0, long n, GEN T, GEN N)
34 {
35 long dP=degpol(P), i, k, m;
36 pari_sp av1, av2;
37 GEN s,y,P_lead;
38
39 if (n<0) pari_err(impl,"polsym of a negative n");
40 if (typ(P) != t_POL) pari_err(typeer,"polsym");
41 if (!signe(P)) pari_err(zeropoler,"polsym");
42 y = cgetg(n+2,t_COL);
43 if (y0)
44 {
45 if (typ(y0) != t_COL) pari_err(typeer,"polsym_gen");
46 m = lg(y0)-1;
47 for (i=1; i<=m; i++) y[i] = y0[i]; /* not memory clean */
48 }
49 else
50 {
51 m = 1;
52 gel(y,1) = stoi(dP);
53 }
54 P += 2; /* strip codewords */
55
56 P_lead = gel(P,dP); if (gcmp1(P_lead)) P_lead = NULL;
57 if (P_lead)
58 {
59 if (N) P_lead = Fq_inv(P_lead,T,N);
60 else if (T) P_lead = QXQ_inv(P_lead,T);
61 }
62 for (k=m; k<=n; k++)
63 {
64 av1 = avma; s = (dP>=k)? gmulsg(k,gel(P,dP-k)): gen_0;
65 for (i=1; i<k && i<=dP; i++)
66 s = gadd(s, gmul(gel(y,k-i+1),gel(P,dP-i)));
67 if (N)
68 {
69 s = Fq_red(s, T, N);
70 if (P_lead) s = Fq_mul(s, P_lead, T, N);
71 }
72 else if (T)
73 {
74 s = grem(s, T);
75 if (P_lead) s = grem(gmul(s, P_lead), T);
76 }
77 else
78 if (P_lead) s = gdiv(s, P_lead);
79 av2 = avma; gel(y,k+1) = gerepile(av1,av2, gneg(s));
80 }
81 return y;
82 }
83
84 GEN
polsym(GEN x,long n)85 polsym(GEN x, long n)
86 {
87 return polsym_gen(x, NULL, n, NULL,NULL);
88 }
89
90 /* assume x and y are polynomials in the same variable whose coeffs can be
91 * compared (used to sort polynomial factorizations)
92 */
93
94 static int
polcmp(void * data,GEN x,GEN y)95 polcmp(void *data, GEN x, GEN y)
96 {
97 int (*coeff_cmp)(GEN,GEN)=(int(*)(GEN,GEN))data;
98 long i, lx = lg(x), ly = lg(y);
99 int fl;
100 if (lx > ly) return 1;
101 if (lx < ly) return -1;
102 for (i=lx-1; i>1; i--)
103 if ((fl = coeff_cmp(gel(x,i), gel(y,i)))) return fl;
104 return 0;
105 }
106
107 /* sort generic factorisation */
108 GEN
sort_factor_gen_aux(GEN y,void * data,int (* cmp)(void *,GEN,GEN))109 sort_factor_gen_aux(GEN y, void *data, int (*cmp)(void *,GEN,GEN))
110 {
111 long n, i;
112 pari_sp av = avma;
113 GEN a,b,A,B,w;
114 a = gel(y,1); n = lg(a); A = new_chunk(n);
115 b = gel(y,2); B = new_chunk(n);
116 w = gen_sort_aux(a, cmp_IND | cmp_C, data, cmp);
117 for (i=1; i<n; i++) { A[i] = a[w[i]]; B[i] = b[w[i]]; }
118 for (i=1; i<n; i++) { a[i] = A[i]; b[i] = B[i]; }
119 avma = av; return y;
120 }
121
122 /* sort generic factorisation */
123 GEN
sort_factor_gen(GEN y,int (* cmp)(GEN,GEN))124 sort_factor_gen(GEN y, int (*cmp)(GEN,GEN))
125 {
126 long n, i;
127 pari_sp av = avma;
128 GEN a,b,A,B,w;
129 a = gel(y,1); n = lg(a); A = new_chunk(n);
130 b = gel(y,2); B = new_chunk(n);
131 w = gen_sort(a, cmp_IND | cmp_C, cmp);
132 for (i=1; i<n; i++) { A[i] = a[w[i]]; B[i] = b[w[i]]; }
133 for (i=1; i<n; i++) { a[i] = A[i]; b[i] = B[i]; }
134 avma = av; return y;
135 }
136
137 GEN
sort_factor(GEN y,int (* cmp)(GEN,GEN))138 sort_factor(GEN y,int (*cmp)(GEN,GEN))
139 {
140 (void)sort_factor_gen_aux(y,(void*)cmp,polcmp);
141 return y;
142 }
143
144 GEN
sort_vecpol_gen(GEN a,int (* cmp)(GEN,GEN))145 sort_vecpol_gen(GEN a, int (*cmp)(GEN,GEN))
146 {
147 long n, i;
148 pari_sp av = avma;
149 GEN A,w;
150 n = lg(a); A = new_chunk(n);
151 w = gen_sort_aux(a, cmp_IND | cmp_C,(void*)cmp, polcmp);
152 for (i=1; i<n; i++) A[i] = a[w[i]];
153 for (i=1; i<n; i++) a[i] = A[i];
154 avma = av; return a;
155 }
156
157 /*In place sort*/
158 GEN
sort_vecpol(GEN y,int (* cmp)(GEN,GEN))159 sort_vecpol(GEN y, int (*cmp)(GEN,GEN))
160 {
161 (void)sort_vecpol_gen(y,cmp);
162 return y;
163 }
164
165 /* centered residue x mod p. po2 = shifti(p, -1) or NULL (euclidean residue) */
166 GEN
centermodii(GEN x,GEN p,GEN po2)167 centermodii(GEN x, GEN p, GEN po2)
168 {
169 GEN y = remii(x, p);
170 switch(signe(y))
171 {
172 case 0: break;
173 case 1: if (po2 && absi_cmp(y,po2) > 0) y = subii(y, p);
174 break;
175 case -1: if (!po2 || absi_cmp(y,po2) > 0) y = addii(y, p);
176 break;
177 }
178 return y;
179 }
180
181 long
s_centermod(long x,ulong pp,ulong pps2)182 s_centermod(long x, ulong pp, ulong pps2)
183 {
184 long y = x % (long)pp;
185 if (y < 0) y += pp;
186 return Fl_center(y, pp,pps2);
187 }
188
189 /* for internal use */
190 GEN
centermod_i(GEN x,GEN p,GEN ps2)191 centermod_i(GEN x, GEN p, GEN ps2)
192 {
193 long i, lx;
194 pari_sp av;
195 GEN y;
196
197 if (!ps2) ps2 = shifti(p,-1);
198 switch(typ(x))
199 {
200 case t_INT: return centermodii(x,p,ps2);
201
202 case t_POL: lx = lg(x);
203 y = cgetg(lx,t_POL); y[1] = x[1];
204 for (i=2; i<lx; i++)
205 {
206 av = avma;
207 gel(y,i) = gerepileuptoint(av, centermodii(gel(x,i),p,ps2));
208 }
209 return normalizepol_i(y, lx);
210
211 case t_COL: lx = lg(x);
212 y = cgetg(lx,t_COL);
213 for (i=1; i<lx; i++) gel(y,i) = centermodii(gel(x,i),p,ps2);
214 return y;
215
216 case t_MAT: lx = lg(x);
217 y = cgetg(lx,t_MAT);
218 for (i=1; i<lx; i++) gel(y,i) = centermod_i(gel(x,i),p,ps2);
219 return y;
220
221 case t_VECSMALL: lx = lg(x);
222 {
223 ulong pp = itou(p), pps2 = itou(ps2);
224 y = cgetg(lx,t_VECSMALL);
225 for (i=1; i<lx; i++) y[i] = s_centermod(x[i], pp, pps2);
226 return y;
227 }
228 }
229 return x;
230 }
231
232 GEN
centermod(GEN x,GEN p)233 centermod(GEN x, GEN p) { return centermod_i(x,p,NULL); }
234
235 /* assume same variables, y normalized, non constant */
236 static GEN
polidivis(GEN x,GEN y,GEN bound)237 polidivis(GEN x, GEN y, GEN bound)
238 {
239 long dx, dy, dz, i, j;
240 pari_sp av;
241 GEN z,p1,y_lead;
242
243 dy=degpol(y);
244 dx=degpol(x);
245 dz=dx-dy; if (dz<0) return NULL;
246 z=cgetg(dz+3,t_POL); z[1] = x[1];
247 x += 2; y += 2; z += 2;
248 y_lead = gel(y,dy);
249 if (gcmp1(y_lead)) y_lead = NULL;
250
251 p1 = gel(x,dx);
252 gel(z,dz) = y_lead? diviiexact(p1,y_lead): icopy(p1);
253 for (i=dx-1; i>=dy; i--)
254 {
255 av = avma; p1 = gel(x,i);
256 for (j=i-dy+1; j<=i && j<=dz; j++)
257 p1 = subii(p1, mulii(gel(z,j),gel(y,i-j)));
258 if (y_lead) p1 = diviiexact(p1,y_lead);
259 if (bound && absi_cmp(p1, bound) > 0) return NULL;
260 p1 = gerepileupto(av,p1);
261 gel(z,i-dy) = p1;
262 }
263 av = avma;
264 for (;; i--)
265 {
266 p1 = gel(x,i);
267 /* we always enter this loop at least once */
268 for (j=0; j<=i && j<=dz; j++)
269 p1 = subii(p1, mulii(gel(z,j),gel(y,i-j)));
270 if (!gcmp0(p1)) return NULL;
271 avma = av;
272 if (!i) break;
273 }
274 return z - 2;
275 }
276
277 /***********************************************************************/
278 /** **/
279 /** QUADRATIC HENSEL LIFT (adapted from V. Shoup's NTL) **/
280 /** **/
281 /***********************************************************************/
282 /* Setup for divide/conquer quadratic Hensel lift
283 * a = set of k t_POL in Z[X] = factors over Fp (T=NULL) or Fp[Y]/(T)
284 * V = set of products of factors built as follows
285 * 1) V[1..k] = initial a
286 * 2) iterate:
287 * append to V the two smallest factors (minimal degree) in a, remove them
288 * from a and replace them by their product [net loss for a = 1 factor]
289 *
290 * W = bezout coeffs W[i]V[i] + W[i+1]V[i+1] = 1
291 *
292 * link[i] = -j if V[i] = a[j]
293 * j if V[i] = V[j] * V[j+1]
294 * Arrays (link, V, W) pre-allocated for 2k - 2 elements */
295 static void
BuildTree(GEN link,GEN V,GEN W,GEN a,GEN T,GEN p)296 BuildTree(GEN link, GEN V, GEN W, GEN a, GEN T, GEN p)
297 {
298 long k = lg(a)-1;
299 long i, j, s, minp, mind;
300
301 for (i=1; i<=k; i++) { V[i] = a[i]; link[i] = -i; }
302
303 for (j=1; j <= 2*k-5; j+=2,i++)
304 {
305 minp = j;
306 mind = degpol(V[j]);
307 for (s=j+1; s<i; s++)
308 if (degpol(V[s]) < mind) { minp = s; mind = degpol(V[s]); }
309
310 lswap(V[j], V[minp]);
311 lswap(link[j], link[minp]);
312
313 minp = j+1;
314 mind = degpol(V[j+1]);
315 for (s=j+2; s<i; s++)
316 if (degpol(V[s]) < mind) { minp = s; mind = degpol(V[s]); }
317
318 lswap(V[j+1], V[minp]);
319 lswap(link[j+1], link[minp]);
320
321 if (T)
322 gel(V,i) = FpXQX_mul(gel(V,j), gel(V,j+1), T, p);
323 else
324 gel(V,i) = FpX_mul(gel(V,j), gel(V,j+1), p);
325 link[i] = j;
326 }
327
328 for (j=1; j <= 2*k-3; j+=2)
329 {
330 GEN d, u, v;
331 if (T)
332 d = FpXQX_extgcd(gel(V,j), gel(V,j+1), T, p, &u, &v);
333 else
334 d = FpX_extgcd(gel(V,j), gel(V,j+1), p, &u, &v);
335 if (degpol(d) > 0) pari_err(talker, "relatively prime polynomials expected");
336 d = gel(d,2);
337 if (!gcmp1(d))
338 {
339 if (typ(d)==t_POL)
340 {
341 d = FpXQ_inv(d, T, p);
342 u = FqX_Fq_mul(u, d, T, p);
343 v = FqX_Fq_mul(v, d, T, p);
344 }
345 else
346 {
347 d = Fp_inv(d, p);
348 u = FpX_Fp_mul(u, d, p);
349 v = FpX_Fp_mul(v, d, p);
350 }
351 }
352 gel(W,j) = u;
353 gel(W,j+1) = v;
354 }
355 }
356
357 /* au + bv = 1 (p0), ab = f (p0). Lift mod p1 = p0 pd (<= p0^2).
358 * If noinv is set, don't lift the inverses u and v */
359 static void
HenselLift(GEN V,GEN W,long j,GEN f,GEN T,GEN pd,GEN p0,int noinv)360 HenselLift(GEN V, GEN W, long j, GEN f, GEN T, GEN pd, GEN p0, int noinv)
361 {
362 pari_sp av = avma;
363 long space = lg(f) * (lgefint(pd) + lgefint(p0));
364 GEN a2,b2,g,z,s,t;
365 GEN a = gel(V,j), b = gel(V,j+1);
366 GEN u = gel(W,j), v = gel(W,j+1);
367
368 if (T) space *= lg(T);
369
370 (void)new_chunk(space); /* HACK */
371 g = gadd(f, gneg_i(gmul(a,b)));
372 if (T) g = FpXQX_red(g, T, mulii(p0,pd));
373 g = gdivexact(g, p0);
374 if (T)
375 {
376 z = FpXQX_mul(v,g, T,pd);
377 t = FpXQX_divrem(z,a, T,pd, &s);
378 }
379 else
380 {
381 g = FpX_red(g, pd);
382 z = FpX_mul(v,g, pd);
383 t = FpX_divrem(z,a, pd, &s);
384 }
385 t = gadd(gmul(u,g), gmul(t,b)); t = T? FpXQX_red(t, T, pd): FpX_red(t, pd);
386 t = gmul(t,p0);
387 s = gmul(s,p0);
388 avma = av;
389
390 /* already reduced mod p1 = pd p0 */
391 a2 = gadd(a,s); gel(V,j) = a2;
392 b2 = gadd(b,t); gel(V,j+1) = b2;
393 if (noinv) return;
394
395 av = avma;
396 (void)new_chunk(space); /* HACK */
397 g = gadd(gmul(u,a2), gmul(v,b2));
398 g = gadd(gneg_i(g), gen_1);
399
400 if (T) g = FpXQX_red(g, T, mulii(p0,pd));
401 g = gdivexact(g, p0);
402 if (T)
403 {
404 z = FpXQX_mul(v,g, T,pd);
405 t = FpXQX_divrem(z,a, T,pd, &s);
406 }
407 else
408 {
409 g = FpX_red(g, pd);
410 z = FpX_mul(v,g, pd);
411 t = FpX_divrem(z,a, pd, &s);
412 }
413 t = gadd(gmul(u,g), gmul(t,b)); t = T? FpXQX_red(t, T, pd): FpX_red(t, pd);
414 t = gmul(t,p0);
415 s = gmul(s,p0);
416 avma = av;
417
418 u = gadd(u,t); gel(W,j) = u;
419 v = gadd(v,s); gel(W,j+1) = v;
420 }
421
422 /* v list of factors, w list of inverses. f = v[j] v[j+1]
423 * Lift v[j] and v[j+1] mod p0 pd (possibly mod T), then all their divisors */
424 static void
RecTreeLift(GEN link,GEN v,GEN w,GEN T,GEN pd,GEN p0,GEN f,long j,int noinv)425 RecTreeLift(GEN link, GEN v, GEN w, GEN T, GEN pd, GEN p0, GEN f, long j, int noinv)
426 {
427 if (j < 0) return;
428
429 HenselLift(v, w, j, f, T, pd, p0, noinv);
430 RecTreeLift(link, v, w, T, pd, p0, gel(v,j) , link[j ], noinv);
431 RecTreeLift(link, v, w, T, pd, p0, gel(v,j+1), link[j+1], noinv);
432 }
433
434 /* lift from p^{e0} to p^{e1} */
435 static void
TreeLift(GEN link,GEN v,GEN w,GEN T,GEN p,long e0,long e1,GEN f,int noinv)436 TreeLift(GEN link, GEN v, GEN w, GEN T, GEN p, long e0, long e1, GEN f, int noinv)
437 {
438 GEN p0 = powiu(p, e0);
439 GEN pd = powiu(p, e1-e0);
440 RecTreeLift(link, v, w, T, pd, p0, f, lgpol(v), noinv);
441 }
442
443 /* Successive accuracies for a quadratic lift.
444 * Eg 9 --> 9,5,3,2,1 instead of 9,8,4,2,1 */
445 GEN
Newton_exponents(long e)446 Newton_exponents(long e)
447 {
448 GEN E = cgetg(BITS_IN_LONG, t_VECSMALL);
449 long l = 1; E[l++] = e;
450 while (e > 1) { e = (e+1)>>1; E[l++] = e; }
451 setlg(E, l); return E;
452 }
453
454 /* a = modular factors of f mod (p,T) [possibly T=NULL], lift to precision p^e0
455 * flag = 0: standard.
456 * flag = 1: return TreeLift structure
457 * flag = 2: a = TreeLift structure, go on lifting, as flag = 1 otherwise */
458 static GEN
MultiLift(GEN f,GEN a,GEN T,GEN p,long e0,long flag)459 MultiLift(GEN f, GEN a, GEN T, GEN p, long e0, long flag)
460 {
461 long l, i, e = e0, k = lg(a) - 1;
462 GEN E, v, w, link;
463 pari_timer Ti;
464
465 if (k < 2 || e < 1) pari_err(talker, "MultiLift: bad args");
466 if (e == 1) return a;
467 if (typ(a[1]) == t_INT) flag = 2;
468 else if (flag == 2) flag = 1;
469
470 E = Newton_exponents(e);
471 e = 1;
472 l = lg(E)-1;
473
474 if (DEBUGLEVEL > 3) TIMERstart(&Ti);
475
476 if (flag != 2)
477 {
478 v = cgetg(2*k - 2 + 1, t_VEC);
479 w = cgetg(2*k - 2 + 1, t_VEC);
480 link=cgetg(2*k - 2 + 1, t_VECSMALL);
481 BuildTree(link, v, w, a, T, p);
482 if (DEBUGLEVEL > 3) msgTIMER(&Ti, "building tree");
483 }
484 else
485 {
486 e = itos(gel(a,1));
487 link = gel(a,2);
488 v = gel(a,3);
489 w = gel(a,4);
490 }
491
492 for (i = l; i > 1; i--) {
493 if (E[i-1] < e) continue;
494 TreeLift(link, v, w, T, p, E[i], E[i-1], f, (flag == 0) && (i == 2));
495 if (DEBUGLEVEL > 3) msgTIMER(&Ti, "lifting to prec %ld", E[i-1]);
496 }
497
498 if (flag)
499 E = mkvec4(stoi(e0), link, v, w);
500 else
501 {
502 E = cgetg(k+1, t_VEC);
503 for (i = 1; i <= 2*k-2; i++)
504 {
505 long t = link[i];
506 if (t < 0) E[-t] = v[i];
507 }
508 }
509 return E;
510 }
511
512 /* Q list of (coprime, monic) factors of pol mod (T,p). Lift mod p^e = pe.
513 * T may be NULL */
514 GEN
hensel_lift_fact(GEN pol,GEN Q,GEN T,GEN p,GEN pe,long e)515 hensel_lift_fact(GEN pol, GEN Q, GEN T, GEN p, GEN pe, long e)
516 {
517 pari_sp av = avma;
518 if (lg(Q) == 2) return mkvec(pol);
519 pol = FqX_normalize(pol, T, pe);
520 return gerepilecopy(av, MultiLift(pol, Q, T, p, e, 0));
521 }
522
523 /* U = NULL treated as 1 */
524 static void
BezoutPropagate(GEN link,GEN v,GEN w,GEN pe,GEN U,GEN f,long j)525 BezoutPropagate(GEN link, GEN v, GEN w, GEN pe, GEN U, GEN f, long j)
526 {
527 GEN Q, R;
528 if (j < 0) return;
529
530 Q = FpX_mul(gel(v,j), gel(w,j), pe);
531 if (U)
532 {
533 Q = FpXQ_mul(Q, U, f, pe);
534 R = FpX_sub(U, Q, pe);
535 }
536 else
537 R = FpX_Fp_add(FpX_neg(Q, pe), gen_1, pe);
538 gel(w,j+1) = Q; /* 0 mod U v[j], 1 mod (1-U) v[j+1] */
539 gel(w,j) = R; /* 1 mod U v[j], 0 mod (1-U) v[j+1] */
540 BezoutPropagate(link, v, w, pe, R, f, link[j ]);
541 BezoutPropagate(link, v, w, pe, Q, f, link[j+1]);
542 }
543
544 /* as above, but return the Bezout coefficients for the lifted modular factors
545 * U[i] = 1 mod Qlift[i]
546 * 0 mod Qlift[j], j != i */
547 GEN
bezout_lift_fact(GEN pol,GEN Q,GEN p,long e)548 bezout_lift_fact(GEN pol, GEN Q, GEN p, long e)
549 {
550 pari_sp av = avma;
551 GEN E, link, v, w, pe;
552 long i, k = lg(Q)-1;
553 if (k == 1) return mkvec(pol);
554 pe = powiu(p, e);
555 pol = FpX_normalize(pol, pe);
556 E = MultiLift(pol, Q, NULL, p, e, 1);
557 link = gel(E,2);
558 v = gel(E,3);
559 w = gel(E,4);
560 BezoutPropagate(link, v, w, pe, NULL, pol, lgpol(v));
561 E = cgetg(k+1, t_VEC);
562 for (i = 1; i <= 2*k-2; i++)
563 {
564 long t = link[i];
565 if (t < 0) E[-t] = w[i];
566 }
567 return gerepilecopy(av, E);
568 }
569
570 /* Front-end for hensel_lift_fact:
571 lift the factorization of pol mod p given by fct to p^exp (if possible) */
572 GEN
polhensellift(GEN pol,GEN fct,GEN p,long exp)573 polhensellift(GEN pol, GEN fct, GEN p, long exp)
574 {
575 GEN p1, p2;
576 long i, j, l;
577 pari_sp av = avma;
578
579 /* we check the arguments */
580 if (typ(pol) != t_POL) pari_err(talker, "not a polynomial in polhensellift");
581 if ((typ(fct) != t_COL && typ(fct) != t_VEC) || (lg(fct) < 3))
582 pari_err(talker, "not a factorization in polhensellift");
583 if (typ(p) != t_INT) pari_err(talker, "not a prime number in polhensellift");
584 if (exp < 1) pari_err(talker, "not a positive exponent in polhensellift");
585
586 l = lg(pol);
587 for (i = 2; i < l; i++)
588 if (typ(pol[i]) != t_INT)
589 pari_err(talker, "not an integral polynomial in polhensellift");
590 p1 = lift(fct); /* make sure the coeffs are integers and not intmods */
591 l = lg(p1);
592 for (i = 1; i < l; i++)
593 {
594 p2 = gel(p1,i);
595 if (typ(p2) != t_POL)
596 {
597 if (typ(p2) != t_INT)
598 pari_err(talker, "not an integral factorization in polhensellift");
599 gel(p1,i) = scalarpol(p2, varn(pol));
600 }
601 }
602
603 /* then we check that pol \equiv \prod f ; f \in fct mod p */
604 p2 = gel(p1,1);
605 for (j = 2; j < l; j++) p2 = FpX_mul(p2, gel(p1,j), p);
606 if (!gcmp0(FpX_sub(pol, p2, p)))
607 pari_err(talker, "not a correct factorization in polhensellift");
608
609 /* finally we check that the elements of fct are coprime mod p */
610 if (!FpX_is_squarefree(pol, p))
611 {
612 for (i = 1; i < l; i++)
613 for (j = 1; j < i; j++)
614 if (degpol(FpX_gcd(gel(p1,i), gel(p1,j), p)))
615 pari_err(talker, "polhensellift: factors %Z and %Z are not coprime",
616 p1[i], p1[j]);
617 }
618 return gerepilecopy(av, hensel_lift_fact(pol,p1,NULL,p,powiu(p,exp),exp));
619 }
620
621 #if 0
622 /* cf Beauzamy et al: upper bound for
623 * lc(x) * [2^(5/8) / pi^(3/8)] e^(1/4n) 2^(n/2) sqrt([x]_2)/ n^(3/8)
624 * where [x]_2 = sqrt(\sum_i=0^n x[i]^2 / binomial(n,i)). One factor has
625 * all coeffs less than then bound */
626 static GEN
627 two_factor_bound(GEN x)
628 {
629 long i, j, n = lg(x) - 3;
630 pari_sp av = avma;
631 GEN *invbin, c, r = cgetr(3), z;
632
633 x += 2; invbin = (GEN*)new_chunk(n+1);
634 z = real_1(3); /* invbin[i] = 1 / binomial(n, i) */
635 for (i=0,j=n; j >= i; i++,j--)
636 {
637 invbin[i] = invbin[j] = z;
638 z = divrs(mulrs(z, i+1), n-i);
639 }
640 z = invbin[0]; /* = 1 */
641 for (i=0; i<=n; i++)
642 {
643 c = gel(x,i); if (!signe(c)) continue;
644 affir(c, r);
645 z = addrr(z, mulrr(gsqr(r), invbin[i]));
646 }
647 z = shiftr(sqrtr(z), n);
648 z = divrr(z, dbltor(pow((double)n, 0.75)));
649 z = grndtoi(sqrtr(z), &i);
650 z = mulii(z, absi(gel(x,n)));
651 return gerepileuptoint(av, shifti(z, 1));
652 }
653 #endif
654
655 /* A | S ==> |a_i| <= binom(d-1, i-1) || S ||_2 + binom(d-1, i) lc(S) */
656 static GEN
Mignotte_bound(GEN S)657 Mignotte_bound(GEN S)
658 {
659 long i, d = degpol(S);
660 GEN lS, C, binlS, bin, N2, p1;
661
662 N2 = sqrtr(QuickNormL2(S,DEFAULTPREC));
663 binlS = bin = vecbinome(d-1);
664 lS = leading_term(S);
665 if (!is_pm1(lS)) binlS = gmul(lS, bin);
666
667 /* i = 0 */
668 C = gel(binlS,1);
669 /* i = d */
670 p1 = N2; if (gcmp(C, p1) < 0) C = p1;
671 for (i = 1; i < d; i++)
672 {
673 p1 = gadd(gmul(gel(bin,i), N2), gel(binlS,i+1));
674 if (gcmp(C, p1) < 0) C = p1;
675 }
676 return C;
677 }
678 /* A | S ==> |a_i|^2 <= 3^{3/2 + d} / (4 \pi d) [P]_2^2,
679 * where [P]_2 is Bombieri's 2-norm */
680 static GEN
Beauzamy_bound(GEN S)681 Beauzamy_bound(GEN S)
682 {
683 const long prec = DEFAULTPREC;
684 long i, d = degpol(S);
685 GEN bin, lS, s, C;
686 bin = vecbinome(d);
687
688 s = real_0(prec);
689 for (i=0; i<=d; i++)
690 {
691 GEN c = gel(S,i+2);
692 if (gcmp0(c)) continue;
693 /* s += P_i^2 / binomial(d,i) */
694 s = addrr(s, gdiv(itor(sqri(c), prec), gel(bin,i+1)));
695 }
696 /* s = [S]_2^2 */
697 C = powrshalf(stor(3,prec), 3 + 2*d); /* 3^{3/2 + d} */
698 C = gdiv(gmul(C, s), gmulsg(4*d, mppi(prec)));
699 lS = absi(leading_term(S));
700 return mulir(lS, sqrtr(C));
701 }
702
703 static GEN
factor_bound(GEN S)704 factor_bound(GEN S)
705 {
706 pari_sp av = avma;
707 GEN a = Mignotte_bound(S);
708 GEN b = Beauzamy_bound(S);
709 if (DEBUGLEVEL>2)
710 {
711 fprintferr("Mignotte bound: %Z\n",a);
712 fprintferr("Beauzamy bound: %Z\n",b);
713 }
714 return gerepileupto(av, ceil_safe(gmin(a, b)));
715 }
716
717 #if 0
718 /* all factors have coeffs less than the bound */
719 static GEN
720 all_factor_bound(GEN x)
721 {
722 long i, n = degpol(x);
723 GEN t, z = gen_0;
724 for (i=2; i<=n+2; i++) z = addii(z, sqri(gel(x,i)));
725 t = absi(gel(x,n+2));
726 z = addii(t, addsi(1, sqrti(z)));
727 z = mulii(z, binomial(stoi(n-1), n>>1));
728 return shifti(mulii(t,z),1);
729 }
730 #endif
731
732 typedef ulong *uGEN;
733
734 /* Naive recombination of modular factors: combine up to maxK modular
735 * factors, degree <= klim and divisible by hint
736 *
737 * target = polynomial we want to factor
738 * famod = array of modular factors. Product should be congruent to
739 * target/lc(target) modulo p^a
740 * For true factors: S1,S2 <= p^b, with b <= a and p^(b-a) < 2^31 */
741 static GEN
cmbf(GEN pol,GEN famod,GEN bound,GEN p,long a,long b,long maxK,long klim,long hint)742 cmbf(GEN pol, GEN famod, GEN bound, GEN p, long a, long b,
743 long maxK, long klim,long hint)
744 {
745 long K = 1, cnt = 1, i,j,k, curdeg, lfamod = lg(famod)-1;
746 ulong spa_b, spa_bs2, Sbound;
747 GEN lc, lcpol, pa = powiu(p,a), pas2 = shifti(pa,-1);
748 uGEN trace1 = (uGEN)cgetg(lfamod+1, t_VECSMALL);
749 uGEN trace2 = (uGEN)cgetg(lfamod+1, t_VECSMALL);
750 GEN ind = cgetg(lfamod+1, t_VECSMALL);
751 GEN degpol = cgetg(lfamod+1, t_VECSMALL);
752 GEN degsofar = cgetg(lfamod+1, t_VECSMALL);
753 GEN listmod = cgetg(lfamod+1, t_COL);
754 GEN fa = cgetg(lfamod+1, t_COL);
755
756 if (maxK < 0) maxK = lfamod-1;
757
758 lc = absi(leading_term(pol));
759 if (is_pm1(lc)) lc = NULL;
760 lcpol = lc? gmul(lc,pol): pol;
761
762 {
763 GEN pa_b,pa_bs2,pb, lc2 = lc? sqri(lc): NULL;
764
765 pa_b = powiu(p, a-b); /* < 2^31 */
766 pa_bs2 = shifti(pa_b,-1);
767 pb= powiu(p, b);
768 for (i=1; i <= lfamod; i++)
769 {
770 GEN T1,T2, P = gel(famod,i);
771 long d = degpol(P);
772
773 degpol[i] = d; P += 2;
774 T1 = gel(P,d-1);/* = - S_1 */
775 T2 = sqri(T1);
776 if (d > 1) T2 = subii(T2, shifti(gel(P,d-2),1));
777 T2 = modii(T2, pa); /* = S_2 Newton sum */
778 if (lc)
779 {
780 T1 = modii(mulii(lc, T1), pa);
781 T2 = modii(mulii(lc2,T2), pa);
782 }
783 trace1[i] = itou(diviiround(T1, pb));
784 trace2[i] = itou(diviiround(T2, pb));
785 }
786 spa_b = (ulong) pa_b[2]; /* < 2^31 */
787 spa_bs2 = (ulong)pa_bs2[2]; /* < 2^31 */
788 }
789 degsofar[0] = 0; /* sentinel */
790
791 /* ind runs through strictly increasing sequences of length K,
792 * 1 <= ind[i] <= lfamod */
793 nextK:
794 if (K > maxK || 2*K > lfamod) goto END;
795 if (DEBUGLEVEL > 3)
796 fprintferr("\n### K = %d, %Z combinations\n", K,binomial(utoipos(lfamod), K));
797 setlg(ind, K+1); ind[1] = 1;
798 Sbound = (ulong) ((K+1)>>1);
799 i = 1; curdeg = degpol[ind[1]];
800 for(;;)
801 { /* try all combinations of K factors */
802 for (j = i; j < K; j++)
803 {
804 degsofar[j] = curdeg;
805 ind[j+1] = ind[j]+1; curdeg += degpol[ind[j+1]];
806 }
807 if (curdeg <= klim && curdeg % hint == 0) /* trial divide */
808 {
809 GEN y, q, list;
810 pari_sp av;
811 ulong t;
812
813 /* d - 1 test */
814 for (t=trace1[ind[1]],i=2; i<=K; i++)
815 t = Fl_add(t, trace1[ind[i]], spa_b);
816 if (t > spa_bs2) t = spa_b - t;
817 if (t > Sbound)
818 {
819 if (DEBUGLEVEL>6) fprintferr(".");
820 goto NEXT;
821 }
822 /* d - 2 test */
823 for (t=trace2[ind[1]],i=2; i<=K; i++)
824 t = Fl_add(t, trace2[ind[i]], spa_b);
825 if (t > spa_bs2) t = spa_b - t;
826 if (t > Sbound)
827 {
828 if (DEBUGLEVEL>6) fprintferr("|");
829 goto NEXT;
830 }
831
832 av = avma;
833 /* check trailing coeff */
834 y = lc;
835 for (i=1; i<=K; i++)
836 {
837 GEN q = constant_term(gel(famod,ind[i]));
838 if (y) q = mulii(y, q);
839 y = centermod_i(q, pa, pas2);
840 }
841 if (!signe(y) || remii(constant_term(lcpol), y) != gen_0)
842 {
843 if (DEBUGLEVEL>3) fprintferr("T");
844 avma = av; goto NEXT;
845 }
846 y = lc; /* full computation */
847 for (i=1; i<=K; i++)
848 {
849 GEN q = gel(famod,ind[i]);
850 if (y) q = gmul(y, q);
851 y = centermod_i(q, pa, pas2);
852 }
853
854 /* y is the candidate factor */
855 if (! (q = polidivis(lcpol,y,bound)) )
856 {
857 if (DEBUGLEVEL>3) fprintferr("*");
858 avma = av; goto NEXT;
859 }
860 /* found a factor */
861 list = cgetg(K+1, t_VEC);
862 gel(listmod,cnt) = list;
863 for (i=1; i<=K; i++) list[i] = famod[ind[i]];
864
865 y = Q_primpart(y);
866 gel(fa,cnt++) = y;
867 /* fix up pol */
868 pol = q;
869 if (lc) pol = Q_div_to_int(pol, leading_term(y));
870 for (i=j=k=1; i <= lfamod; i++)
871 { /* remove used factors */
872 if (j <= K && i == ind[j]) j++;
873 else
874 {
875 famod[k] = famod[i];
876 trace1[k] = trace1[i];
877 trace2[k] = trace2[i];
878 degpol[k] = degpol[i]; k++;
879 }
880 }
881 lfamod -= K;
882 if (lfamod < 2*K) goto END;
883 i = 1; curdeg = degpol[ind[1]];
884 bound = factor_bound(pol);
885 if (lc) lc = absi(leading_term(pol));
886 lcpol = lc? gmul(lc,pol): pol;
887 if (DEBUGLEVEL>3)
888 fprintferr("\nfound factor %Z\nremaining modular factor(s): %ld\n",
889 y, lfamod);
890 continue;
891 }
892
893 NEXT:
894 for (i = K+1;;)
895 {
896 if (--i == 0) { K++; goto nextK; }
897 if (++ind[i] <= lfamod - K + i)
898 {
899 curdeg = degsofar[i-1] + degpol[ind[i]];
900 if (curdeg <= klim) break;
901 }
902 }
903 }
904 END:
905 if (degpol(pol) > 0)
906 { /* leftover factor */
907 if (signe(leading_term(pol)) < 0) pol = gneg_i(pol);
908
909 setlg(famod, lfamod+1);
910 gel(listmod,cnt) = shallowcopy(famod);
911 gel(fa,cnt++) = pol;
912 }
913 if (DEBUGLEVEL>6) fprintferr("\n");
914 setlg(listmod, cnt);
915 setlg(fa, cnt); return mkvec2(fa, listmod);
916 }
917
918 void
factor_quad(GEN x,GEN res,long * ptcnt)919 factor_quad(GEN x, GEN res, long *ptcnt)
920 {
921 GEN a = gel(x,4), b = gel(x,3), c = gel(x,2), d, u, z1, z2, t;
922 GEN D = subii(sqri(b), shifti(mulii(a,c), 2));
923 long v, cnt = *ptcnt;
924
925 if (!Z_issquarerem(D, &d)) { gel(res,cnt++) = x; *ptcnt = cnt; return; }
926
927 t = shifti(negi(addii(b, d)), -1);
928 z1 = gdiv(t, a); u = denom(z1);
929 z2 = gdiv(addii(t, d), a);
930 v = varn(x);
931 gel(res,cnt++) = gmul(u, gsub(pol_x[v], z1)); u = diviiexact(a, u);
932 gel(res,cnt++) = gmul(u, gsub(pol_x[v], z2)); *ptcnt = cnt;
933 }
934
935 /* y > 1 and B integers. Let n such that y^(n-1) <= B < y^n.
936 * Return e = max(n,1), set *ptq = y^e if non-NULL */
937 long
logint(GEN B,GEN y,GEN * ptq)938 logint(GEN B, GEN y, GEN *ptq)
939 {
940 pari_sp av = avma;
941 long e,i,fl;
942 GEN q,pow2, r = y;
943
944 if (typ(B) != t_INT) B = ceil_safe(B);
945 if (expi(B) <= (expi(y) << 6)) /* e small, be naive */
946 {
947 for (e=1; cmpii(r, B) < 0; e++) r = mulii(r,y);
948 goto END;
949 }
950 /* binary splitting: compute bits of e one by one */
951 /* compute pow2[i] = y^(2^i) [i < very crude upper bound for log_2(n)] */
952 pow2 = new_chunk(bit_accuracy(lgefint(B)));
953 gel(pow2,0) = y;
954 for (i=0,q=r;; )
955 {
956 fl = cmpii(r,B); if (fl >= 0) break;
957 q = r; r = sqri(q);
958 i++; gel(pow2,i) = r;
959 }
960 if (i == 0) { e = 1; goto END; } /* y <= B */
961
962 for (i--, e=1<<i;;)
963 { /* y^e = q < B <= r = q * y^(2^i) */
964 if (!fl) break; /* B = r */
965 /* q < B < r */
966 if (--i < 0) { if (fl > 0) e++; break; }
967 r = mulii(q, gel(pow2,i));
968 fl = cmpii(r, B);
969 if (fl <= 0) { e += (1<<i); q = r; }
970 }
971 if (fl <= 0) { e++; r = mulii(r,y); }
972 END:
973 if (ptq) *ptq = gerepileuptoint(av, icopy(r)); else avma = av;
974 return e;
975 }
976
977 /* recombination of modular factors: van Hoeij's algorithm */
978
979 /* Q in Z[X], return Q(2^n) */
980 static GEN
shifteval(GEN Q,long n)981 shifteval(GEN Q, long n)
982 {
983 long i, l = lg(Q);
984 GEN s;
985
986 if (!signe(Q)) return gen_0;
987 s = gel(Q,l-1);
988 for (i = l-2; i > 1; i--) s = addii(gel(Q,i), shifti(s, n));
989 return s;
990 }
991
992 /* return integer y such that all |a| <= y if P(a) = 0 */
993 static GEN
root_bound(GEN P0)994 root_bound(GEN P0)
995 {
996 GEN Q = shallowcopy(P0), lP = absi(leading_term(Q)), x,y,z;
997 long k, d = degpol(Q);
998
999 /* P0 = lP x^d + Q, deg Q < d */
1000 Q = normalizepol_i(Q, d+2);
1001 for (k=lg(Q)-1; k>1; k--) gel(Q,k) = absi(gel(Q,k));
1002 k = (long)(cauchy_bound(P0) / LOG2);
1003 for ( ; k >= 0; k--)
1004 {
1005 pari_sp av = avma;
1006 /* y = 2^k; Q(y) >= lP y^d ? */
1007 if (cmpii(shifteval(Q,k), shifti(lP, d*k)) >= 0) break;
1008 avma = av;
1009 }
1010 if (k < 0) k = 0;
1011 x = int2n(k);
1012 y = int2n(k+1);
1013 for(k=0; ; k++)
1014 {
1015 z = shifti(addii(x,y), -1);
1016 if (equalii(x,z) || k > 5) break;
1017 if (cmpii(poleval(Q,z), mulii(lP, powiu(z, d))) < 0)
1018 y = z;
1019 else
1020 x = z;
1021 }
1022 return y;
1023 }
1024
1025 static GEN
ZM_HNFimage(GEN x)1026 ZM_HNFimage(GEN x)
1027 {
1028 return (lg(x) > 50)? hnflll_i(x, NULL, 1): hnfall_i(x, NULL, 1);
1029 }
1030
1031 GEN
special_pivot(GEN x)1032 special_pivot(GEN x)
1033 {
1034 GEN t, H = ZM_HNFimage(x);
1035 long i,j, l = lg(H), h = lg(H[1]);
1036 for (i=1; i<h; i++)
1037 {
1038 int fl = 0;
1039 for (j=1; j<l; j++)
1040 {
1041 t = gcoeff(H,i,j);
1042 if (signe(t))
1043 {
1044 if (!is_pm1(t) || fl) return NULL;
1045 fl = 1;
1046 }
1047 }
1048 }
1049 return H;
1050 }
1051
1052 /* B from lllint_i: return [ |b_i^*|^2, i ] */
1053 GEN
GS_norms(GEN B,long prec)1054 GS_norms(GEN B, long prec)
1055 {
1056 long i, l = lg(B);
1057 GEN v = gmul(B, real_1(prec));
1058 l--; setlg(v, l);
1059 for (i=1; i<l; i++)
1060 gel(v,i) = divrr(gel(v,i+1), gel(v,i));
1061 return v;
1062 }
1063
1064 GEN
chk_factors_get(GEN lt,GEN famod,GEN c,GEN T,GEN N)1065 chk_factors_get(GEN lt, GEN famod, GEN c, GEN T, GEN N)
1066 {
1067 long i = 1, j, l = lg(famod);
1068 GEN V = cgetg(l, t_VEC);
1069 for (j = 1; j < l; j++)
1070 if (signe(c[j])) V[i++] = famod[j];
1071 if (lt && i > 1) gel(V,1) = gmul(lt, gel(V,1));
1072 setlg(V, i);
1073 if (T) return FpXQXV_prod(V, T, N);
1074 else return FpXV_prod(V,N);
1075 }
1076
1077 static GEN
chk_factors(GEN P,GEN M_L,GEN bound,GEN famod,GEN pa)1078 chk_factors(GEN P, GEN M_L, GEN bound, GEN famod, GEN pa)
1079 {
1080 long i, r;
1081 GEN pol = P, list, piv, y, ltpol, lt;
1082
1083 piv = special_pivot(M_L);
1084 if (!piv) return NULL;
1085 if (DEBUGLEVEL>3) fprintferr("special_pivot output:\n%Z\n",piv);
1086
1087 r = lg(piv)-1;
1088 list = cgetg(r+1, t_COL);
1089 lt = absi(leading_term(pol));
1090 if (is_pm1(lt)) lt = NULL;
1091 ltpol = lt? gmul(lt, pol): pol;
1092 for (i = 1;;)
1093 {
1094 if (DEBUGLEVEL) fprintferr("LLL_cmbf: checking factor %ld\n",i);
1095 y = chk_factors_get(lt, famod, gel(piv,i), NULL, pa);
1096 y = FpX_center(y, pa);
1097 if (! (pol = polidivis(ltpol,y,bound)) ) return NULL;
1098 if (lt) y = Q_primpart(y);
1099 gel(list,i) = y;
1100 if (++i >= r) break;
1101
1102 if (lt)
1103 {
1104 pol = gdivexact(pol, leading_term(y));
1105 lt = absi(leading_term(pol));
1106 ltpol = gmul(lt, pol);
1107 }
1108 else
1109 ltpol = pol;
1110 }
1111 y = Q_primpart(pol);
1112 gel(list,i) = y; return list;
1113 }
1114
1115 GEN
LLL_check_progress(GEN Bnorm,long n0,GEN m,int final,long * ti_LLL)1116 LLL_check_progress(GEN Bnorm, long n0, GEN m, int final, long *ti_LLL)
1117 {
1118 GEN B, norm, u;
1119 long i, R;
1120 pari_timer T;
1121
1122 if (DEBUGLEVEL>2) TIMERstart(&T);
1123 u = lllint_i(m, final? 1000: 4, 0, NULL, NULL, &B);
1124 if (DEBUGLEVEL>2) *ti_LLL += TIMER(&T);
1125 norm = GS_norms(B, DEFAULTPREC);
1126 for (R=lg(m)-1; R > 0; R--)
1127 if (cmprr(gel(norm,R), Bnorm) < 0) break;
1128 for (i=1; i<=R; i++) setlg(u[i], n0+1);
1129 if (R <= 1)
1130 {
1131 if (!R) pari_err(bugparier,"LLL_cmbf [no factor]");
1132 return NULL; /* irreducible */
1133 }
1134 setlg(u, R+1); return u;
1135 }
1136
1137 static ulong
next2pow(ulong a)1138 next2pow(ulong a)
1139 {
1140 ulong b = 1;
1141 while (b < a) b <<= 1;
1142 return b;
1143 }
1144
1145 /* Recombination phase of Berlekamp-Zassenhaus algorithm using a variant of
1146 * van Hoeij's knapsack
1147 *
1148 * P = squarefree in Z[X].
1149 * famod = array of (lifted) modular factors mod p^a
1150 * bound = Mignotte bound for the size of divisors of P (for the sup norm)
1151 * previously recombined all set of factors with less than rec elts */
1152 static GEN
LLL_cmbf(GEN P,GEN famod,GEN p,GEN pa,GEN bound,long a,long rec)1153 LLL_cmbf(GEN P, GEN famod, GEN p, GEN pa, GEN bound, long a, long rec)
1154 {
1155 const long N0 = 1; /* # of traces added at each step */
1156 double BitPerFactor = 0.5; /* nb bits in p^(a-b) / modular factor */
1157 long i,j,tmax,n0,C, dP = degpol(P);
1158 double logp = log((double)itos(p)), LOGp2 = LOG2/logp;
1159 double b0 = log((double)dP*2) / logp, logBr;
1160 GEN lP, Br, Bnorm, Tra, T2, TT, CM_L, m, list, ZERO;
1161 pari_sp av, av2, lim;
1162 long ti_LLL = 0, ti_CF = 0;
1163
1164 lP = absi(leading_term(P));
1165 if (is_pm1(lP)) lP = NULL;
1166 Br = root_bound(P);
1167 if (lP) Br = gmul(lP, Br);
1168 logBr = gtodouble(glog(Br, DEFAULTPREC)) / logp;
1169
1170 n0 = lg(famod) - 1;
1171 C = (long)ceil( sqrt(N0 * n0 / 4.) ); /* > 1 */
1172 Bnorm = dbltor(n0 * (C*C + N0*n0/4.) * 1.00001);
1173 ZERO = zeromat(n0, N0);
1174
1175 av = avma; lim = stack_lim(av, 1);
1176 TT = cgetg(n0+1, t_VEC);
1177 Tra = cgetg(n0+1, t_MAT);
1178 for (i=1; i<=n0; i++)
1179 {
1180 TT[i] = 0;
1181 gel(Tra,i) = cgetg(N0+1, t_COL);
1182 }
1183 CM_L = gscalsmat(C, n0);
1184 /* tmax = current number of traces used (and computed so far) */
1185 for (tmax = 0;; tmax += N0)
1186 {
1187 long b, bmin, bgood, delta, tnew = tmax + N0, r = lg(CM_L)-1;
1188 GEN M_L, q, CM_Lp, oldCM_L;
1189 int first = 1;
1190 pari_timer ti2, TI;
1191
1192 bmin = (long)ceil(b0 + tnew*logBr);
1193 if (DEBUGLEVEL>2)
1194 fprintferr("\nLLL_cmbf: %ld potential factors (tmax = %ld, bmin = %ld)\n",
1195 r, tmax, bmin);
1196
1197 /* compute Newton sums (possibly relifting first) */
1198 if (a <= bmin)
1199 {
1200 a = (long)ceil(bmin + 3*N0*logBr) + 1; /* enough for 3 more rounds */
1201 a = (long)next2pow((ulong)a);
1202
1203 pa = powiu(p,a);
1204 famod = hensel_lift_fact(P,famod,NULL,p,pa,a);
1205 for (i=1; i<=n0; i++) TT[i] = 0;
1206 }
1207 for (i=1; i<=n0; i++)
1208 {
1209 GEN p1 = gel(Tra,i);
1210 GEN p2 = polsym_gen(gel(famod,i), gel(TT,i), tnew, NULL, pa);
1211 gel(TT,i) = p2;
1212 p2 += 1+tmax; /* ignore traces number 0...tmax */
1213 for (j=1; j<=N0; j++) p1[j] = p2[j];
1214 if (lP)
1215 { /* make Newton sums integral */
1216 GEN lPpow = powiu(lP, tmax);
1217 for (j=1; j<=N0; j++)
1218 {
1219 lPpow = mulii(lPpow,lP);
1220 gel(p1,j) = mulii(gel(p1,j), lPpow);
1221 }
1222 }
1223 }
1224
1225 /* compute truncation parameter */
1226 if (DEBUGLEVEL>2) { TIMERstart(&ti2); TIMERstart(&TI); }
1227 oldCM_L = CM_L;
1228 av2 = avma;
1229 delta = b = 0; /* -Wall */
1230 AGAIN:
1231 M_L = Q_div_to_int(CM_L, utoipos(C));
1232 T2 = centermod( gmul(Tra, M_L), pa );
1233 if (first)
1234 { /* initialize lattice, using few p-adic digits for traces */
1235 double t = gexpo(T2) - max(32, BitPerFactor*r);
1236 bgood = (long) (t * LOGp2);
1237 b = max(bmin, bgood);
1238 delta = a - b;
1239 }
1240 else
1241 { /* add more p-adic digits and continue reduction */
1242 long b0 = (long)(gexpo(T2) * LOGp2);
1243 if (b0 < b) b = b0;
1244 b = max(b-delta, bmin);
1245 if (b - delta/2 < bmin) b = bmin; /* near there. Go all the way */
1246 }
1247
1248 q = powiu(p, b);
1249 m = vconcat( CM_L, gdivround(T2, q) );
1250 if (first)
1251 {
1252 GEN P1 = gscalmat(powiu(p, a-b), N0);
1253 first = 0;
1254 m = shallowconcat( m, vconcat(ZERO, P1) );
1255 /* [ C M_L 0 ]
1256 * m = [ ] square matrix
1257 * [ T2' p^(a-b) I_N0 ] T2' = Tra * M_L truncated
1258 */
1259 }
1260
1261 CM_L = LLL_check_progress(Bnorm, n0, m, b == bmin, /*dbg:*/ &ti_LLL);
1262 if (DEBUGLEVEL>2)
1263 fprintferr("LLL_cmbf: (a,b) =%4ld,%4ld; r =%3ld -->%3ld, time = %ld\n",
1264 a,b, lg(m)-1, CM_L? lg(CM_L)-1: 1, TIMER(&TI));
1265 if (!CM_L) { list = mkcol(P); break; }
1266 if (b > bmin)
1267 {
1268 CM_L = gerepilecopy(av2, CM_L);
1269 goto AGAIN;
1270 }
1271 if (DEBUGLEVEL>2) msgTIMER(&ti2, "for this block of traces");
1272
1273 i = lg(CM_L) - 1;
1274 if (i == r && gequal(CM_L, oldCM_L))
1275 {
1276 CM_L = oldCM_L;
1277 avma = av2; continue;
1278 }
1279
1280 CM_Lp = FpM_image(CM_L, utoipos(27449)); /* inexpensive test */
1281 if (lg(CM_Lp) != lg(CM_L))
1282 {
1283 if (DEBUGLEVEL>2) fprintferr("LLL_cmbf: rank decrease\n");
1284 CM_L = ZM_HNFimage(CM_L);
1285 }
1286
1287 if (i <= r && i*rec < n0)
1288 {
1289 pari_timer ti;
1290 if (DEBUGLEVEL>2) TIMERstart(&ti);
1291 list = chk_factors(P, Q_div_to_int(CM_L,utoipos(C)), bound, famod, pa);
1292 if (DEBUGLEVEL>2) ti_CF += TIMER(&ti);
1293 if (list) break;
1294 if (DEBUGLEVEL>2) fprintferr("LLL_cmbf: chk_factors failed");
1295 CM_L = gerepilecopy(av2, CM_L);
1296 }
1297 if (low_stack(lim, stack_lim(av,1)))
1298 {
1299 if(DEBUGMEM>1) pari_warn(warnmem,"LLL_cmbf");
1300 gerepileall(av, 5, &CM_L, &TT, &Tra, &famod, &pa);
1301 }
1302 }
1303 if (DEBUGLEVEL>2)
1304 fprintferr("* Time LLL: %ld\n* Time Check Factor: %ld\n",ti_LLL,ti_CF);
1305 return list;
1306 }
1307
1308 /* Find a,b minimal such that A < q^a, B < q^b, 1 << q^(a-b) < 2^31 */
1309 int
cmbf_precs(GEN q,GEN A,GEN B,long * pta,long * ptb,GEN * qa,GEN * qb)1310 cmbf_precs(GEN q, GEN A, GEN B, long *pta, long *ptb, GEN *qa, GEN *qb)
1311 {
1312 long a,b,amin,d = (long)(31 * LOG2/gtodouble(glog(q,DEFAULTPREC)) - 1e-5);
1313 int fl = 0;
1314
1315 b = logint(B, q, qb);
1316 amin = b + d;
1317 if (gcmp(powiu(q, amin), A) <= 0)
1318 {
1319 a = logint(A, q, qa);
1320 b = a - d; *qb = powiu(q, b);
1321 }
1322 else
1323 { /* not enough room */
1324 a = amin; *qa = powiu(q, a);
1325 fl = 1;
1326 }
1327 if (DEBUGLEVEL > 3) {
1328 fprintferr("S_2 bound: %Z^%ld\n", q,b);
1329 fprintferr("coeff bound: %Z^%ld\n", q,a);
1330 }
1331 *pta = a;
1332 *ptb = b; return fl;
1333 }
1334
1335 /* use van Hoeij's knapsack algorithm */
1336 GEN
combine_factors(GEN target,GEN famod,GEN p,long klim,long hint)1337 combine_factors(GEN target, GEN famod, GEN p, long klim, long hint)
1338 {
1339 GEN la, B, A, res, L, pa, pb, listmod;
1340 long a,b, l, maxK = 3, nft = lg(famod)-1, n = degpol(target);
1341 pari_timer T;
1342
1343 A = factor_bound(target);
1344
1345 la = absi(leading_term(target));
1346 B = mulsi(n, sqri(gmul(la, root_bound(target)))); /* = bound for S_2 */
1347
1348 (void)cmbf_precs(p, A, B, &a, &b, &pa, &pb);
1349
1350 if (DEBUGLEVEL>2) (void)TIMER(&T);
1351 famod = hensel_lift_fact(target,famod,NULL,p,pa,a);
1352 if (nft < 11) maxK = -1; /* few modular factors: try all posibilities */
1353 if (DEBUGLEVEL>2) msgTIMER(&T, "Hensel lift (mod %Z^%ld)", p,a);
1354 L = cmbf(target, famod, A, p, a, b, maxK, klim, hint);
1355 if (DEBUGLEVEL>2) msgTIMER(&T, "Naive recombination");
1356
1357 res = gel(L,1);
1358 listmod = gel(L,2); l = lg(listmod)-1;
1359 famod = gel(listmod,l);
1360 if (maxK >= 0 && lg(famod)-1 > 2*maxK)
1361 {
1362 if (l!=1) A = factor_bound(gel(res,l));
1363 if (DEBUGLEVEL > 4) fprintferr("last factor still to be checked\n");
1364 L = LLL_cmbf(gel(res,l), famod, p, pa, A, a, maxK);
1365 if (DEBUGLEVEL>2) msgTIMER(&T,"Knapsack");
1366 /* remove last elt, possibly unfactored. Add all new ones. */
1367 setlg(res, l); res = shallowconcat(res, L);
1368 }
1369 return res;
1370 }
1371
1372 /* assume pol(0) != 0, polp = pol/lc(pol) mod p.
1373 * Return vector of rational roots of a */
1374 GEN
DDF_roots(GEN pol,GEN polp,GEN p)1375 DDF_roots(GEN pol, GEN polp, GEN p)
1376 {
1377 GEN lc, lcpol, z, pe, pes2, bound;
1378 long i, m, e, lz, v = varn(pol);
1379 pari_sp av, lim;
1380 pari_timer T;
1381
1382 if (DEBUGLEVEL>2) TIMERstart(&T);
1383 lc = absi(leading_term(pol));
1384 if (is_pm1(lc)) lc = NULL;
1385 lcpol = lc? gmul(lc,pol): pol;
1386
1387 bound = root_bound(pol);
1388 if (lc) bound = mulii(lc, bound);
1389 e = logint(addis(shifti(bound, 1), 1), p, &pe);
1390 pes2 = shifti(pe, -1);
1391 if (DEBUGLEVEL>2) msgTIMER(&T, "Root bound");
1392
1393 av = avma; lim = stack_lim(av,2);
1394 z = FpX_roots(polp, p);
1395 lz = lg(z)-1;
1396 if (lz > (degpol(pol) >> 2))
1397 { /* many roots */
1398 z = shallowconcat(deg1_from_roots(z, v),
1399 FpX_div(polp, FpV_roots_to_pol(z, p, v), p));
1400 z = hensel_lift_fact(pol, z, NULL, p, pe, e);
1401 }
1402 else
1403 {
1404 z = ZpX_liftroots(pol, z, p, e);
1405 z = deg1_from_roots(z, v);
1406 }
1407 if (DEBUGLEVEL>2) msgTIMER(&T, "Hensel lift (mod %Z^%ld)", p,e);
1408
1409 for (m=1, i=1; i <= lz; i++)
1410 {
1411 GEN q, r, y = gel(z,i);
1412 if (lc) y = gmul(y, lc);
1413 y = centermod_i(y, pe, pes2);
1414 if (! (q = polidivis(lcpol, y, NULL)) ) continue;
1415
1416 lcpol = pol = q;
1417 r = negi( constant_term(y) );
1418 if (lc) {
1419 r = gdiv(r,lc);
1420 pol = Q_primpart(pol);
1421 lc = absi( leading_term(pol) );
1422 if (is_pm1(lc)) lc = NULL; else lcpol = gmul(lc, pol);
1423 }
1424 gel(z,m++) = r;
1425 if (low_stack(lim, stack_lim(av,2)))
1426 {
1427 if (DEBUGMEM>1) pari_warn(warnmem,"DDF_roots, m = %ld", m);
1428 gerepileall(av, lc? 4:2, &z, &pol, &lc, &lcpol);
1429
1430 }
1431 }
1432 if (DEBUGLEVEL>2) msgTIMER(&T, "Recombination");
1433 z[0] = evaltyp(t_VEC) | evallg(m); return z;
1434 }
1435
1436 /* Assume a squarefree, degree(a) > 0, a(0) != 0.
1437 * If fl != 0 look only for rational roots */
1438 static GEN
DDF(GEN a,long hint,int fl)1439 DDF(GEN a, long hint, int fl)
1440 {
1441 GEN lead, prime, famod, z, ap;
1442 const long da = degpol(a);
1443 long chosenp, p, nfacp, np, nmax, ti = 0;
1444 pari_sp av = avma, av1;
1445 byteptr pt = diffptr;
1446 const long MAXNP = 7;
1447 pari_timer T, T2;
1448
1449 if (DEBUGLEVEL>2) { TIMERstart(&T); TIMERstart(&T2); }
1450 if (hint <= 0) hint = 1;
1451 nmax = da+1;
1452 chosenp = 0;
1453 lead = gel(a,da+2); if (gcmp1(lead)) lead = NULL;
1454 av1 = avma;
1455 for (p = np = 0; np < MAXNP; avma = av1)
1456 {
1457 NEXT_PRIME_VIADIFF_CHECK(p,pt);
1458 if (lead && !smodis(lead,p)) continue;
1459 z = ZX_to_Flx(a, p);
1460 if (!Flx_is_squarefree(z, p)) continue;
1461
1462 nfacp = fl? Flx_nbroots(z, p): Flx_nbfact(z, p);
1463 if (DEBUGLEVEL>4)
1464 fprintferr("...tried prime %3ld (%-3ld %s). Time = %ld\n",
1465 p, nfacp, fl?"roots": "factors", TIMER(&T2));
1466 if (nfacp < nmax)
1467 {
1468 if (nfacp <= 1)
1469 {
1470 if (!fl) { avma = av; return mkcol(a); } /* irreducible */
1471 if (!nfacp) return cgetg(1, t_VEC); /* no root */
1472 }
1473 nmax = nfacp; chosenp = p;
1474 if (da > 100 && nmax < 5) break; /* large degree, few factors. Enough */
1475 }
1476 np++;
1477 }
1478 prime = utoipos(chosenp);
1479 ap = lead? FpX_normalize(a, prime): FpX_red(a, prime);
1480 if (fl) return gerepilecopy(av, DDF_roots(a, ap, prime));
1481
1482 famod = cgetg(nmax+1,t_COL);
1483 gel(famod,1) = ap;
1484 if (nmax != FpX_split_Berlekamp((GEN*)(famod+1), prime))
1485 pari_err(bugparier,"DDF: wrong numbers of factors");
1486 if (DEBUGLEVEL>2)
1487 {
1488 if (DEBUGLEVEL>4) msgTIMER(&T2, "splitting mod p = %ld", chosenp);
1489 ti = TIMER(&T);
1490 fprintferr("Time setup: %ld\n", ti);
1491 }
1492 z = combine_factors(a, famod, prime, da-1, hint);
1493 if (DEBUGLEVEL>2)
1494 fprintferr("Total Time: %ld\n===========\n", ti + TIMER(&T));
1495 return gerepilecopy(av, z);
1496 }
1497
1498 /* A(X^d) --> A(X) */
1499 GEN
poldeflate_i(GEN x0,long d)1500 poldeflate_i(GEN x0, long d)
1501 {
1502 GEN z, y, x;
1503 long i,id, dy, dx = degpol(x0);
1504 if (d <= 1) return x0;
1505 if (dx < 0) return zeropol(varn(x0));
1506 dy = dx/d;
1507 y = cgetg(dy+3, t_POL); y[1] = x0[1];
1508 z = y + 2;
1509 x = x0+ 2;
1510 for (i=id=0; i<=dy; i++,id+=d) z[i] = x[id];
1511 return y;
1512 }
1513
1514 long
checkdeflate(GEN x)1515 checkdeflate(GEN x)
1516 {
1517 long d = 0, i, lx = lg(x);
1518 for (i=3; i<lx; i++)
1519 if (!gcmp0(gel(x,i))) { d = cgcd(d,i-2); if (d == 1) break; }
1520 return d;
1521 }
1522
1523 GEN
gdeflate(GEN x,long v,long d)1524 gdeflate(GEN x, long v, long d)
1525 {
1526 long i, lx, tx = typ(x);
1527 GEN z;
1528 if (is_scalar_t(tx)) return gcopy(x);
1529 if (d <= 0) pari_err(talker,"need positive degree in gdeflate");
1530 if (tx == t_POL || tx == t_SER)
1531 {
1532 long vx = varn(x);
1533 pari_sp av;
1534 if (varncmp(vx, v) < 0)
1535 {
1536 lx = lg(x); z = cgetg(lx, tx); z[1] = x[1];
1537 for (i=2; i<lx; i++) gel(z,i) = gdeflate(gel(x,i),v,d);
1538 return z;
1539 }
1540 if (varncmp(vx, v) > 0) return gcopy(x);
1541 av = avma;
1542 if (tx == t_SER)
1543 {
1544 long V = valp(x);
1545 GEN y;
1546
1547 lx = lg(x);
1548 if (lx == 2) return zeroser(v, V / d);
1549 y = ser2pol_i(x, lx);
1550 if (V % d != 0 || checkdeflate(y) % d != 0)
1551 pari_err(talker, "can't deflate this power series (d = %ld): %Z", d, x);
1552 y = poltoser(poldeflate_i(y, d), v, 1 + (lx-3)/d);
1553 setvalp(y, V/d); return gerepilecopy(av, y);
1554 }
1555 if (checkdeflate(x) % d != 0) pari_err(cant_deflate);
1556 return gerepilecopy(av, poldeflate_i(x,d));
1557 }
1558 if (tx == t_RFRAC)
1559 {
1560 z = cgetg(3, t_RFRAC);
1561 gel(z,1) = gdeflate(gel(x,1),v,d);
1562 gel(z,2) = gdeflate(gel(x,2),v,d);
1563 return z;
1564 }
1565 if (is_matvec_t(tx))
1566 {
1567 lx = lg(x); z = cgetg(lx, tx);
1568 for (i=1; i<lx; i++) gel(z,i) = gdeflate(gel(x,i),v,d);
1569 return z;
1570 }
1571 pari_err(typeer,"gdeflate");
1572 return NULL; /* not reached */
1573 }
1574
1575 /* set *m to the largest d such that x0 = A(X^d); return A */
1576 GEN
poldeflate(GEN x,long * m)1577 poldeflate(GEN x, long *m)
1578 {
1579 *m = checkdeflate(x);
1580 return poldeflate_i(x, *m);
1581 }
1582
1583 /* return x0(X^d) */
1584 GEN
polinflate(GEN x0,long d)1585 polinflate(GEN x0, long d)
1586 {
1587 long i, id, dy, dx = degpol(x0);
1588 GEN x = x0 + 2, z, y;
1589 dy = dx*d;
1590 y = cgetg(dy+3, t_POL); y[1] = x0[1];
1591 z = y + 2;
1592 for (i=0; i<=dy; i++) gel(z,i) = gen_0;
1593 for (i=id=0; i<=dx; i++,id+=d) z[id] = x[i];
1594 return y;
1595 }
1596
1597 /* Distinct Degree Factorization (deflating first)
1598 * Assume x squarefree, degree(x) > 0, x(0) != 0 */
1599 GEN
ZX_DDF(GEN x,long hint)1600 ZX_DDF(GEN x, long hint)
1601 {
1602 GEN L;
1603 long m;
1604 x = poldeflate(x, &m);
1605 L = DDF(x, hint, 0);
1606 if (m > 1)
1607 {
1608 GEN e, v, fa = factoru(m);
1609 long i,j,k, l;
1610
1611 e = gel(fa,2); k = 0;
1612 fa= gel(fa,1); l = lg(fa);
1613 for (i=1; i<l; i++) k += e[i];
1614 v = cgetg(k+1, t_VECSMALL); k = 1;
1615 for (i=1; i<l; i++)
1616 for (j=1; j<=e[i]; j++) v[k++] = fa[i];
1617 for (k--; k; k--)
1618 {
1619 GEN L2 = cgetg(1,t_VEC);
1620 for (i=1; i < lg(L); i++)
1621 L2 = shallowconcat(L2, DDF(polinflate(gel(L,i), v[k]), hint, 0));
1622 L = L2;
1623 }
1624 }
1625 return L;
1626 }
1627
1628 /* SquareFree Factorization. f = prod P^e, all e distinct, in Z[X] (char 0
1629 * would be enough, if modulargcd --> ggcd). Return (P), set *ex = (e) */
1630 GEN
ZX_squff(GEN f,GEN * ex)1631 ZX_squff(GEN f, GEN *ex)
1632 {
1633 GEN T, V, W, P, e;
1634 long i, k, dW, n, val;
1635
1636 if (signe(leading_term(f)) < 0) f = gneg_i(f);
1637 val = ZX_valuation(f, &f);
1638 n = 1 + degpol(f); if (val) n++;
1639 e = cgetg(n,t_VECSMALL);
1640 P = cgetg(n,t_COL);
1641
1642 T = modulargcd(derivpol(f), f);
1643 V = RgX_div(f,T);
1644 for (k=i=1;; k++)
1645 {
1646 W = modulargcd(T,V); T = RgX_div(T,W); dW = degpol(W);
1647 /* W = prod P^e, e > k; V = prod P^e, e >= k */
1648 if (dW != degpol(V)) { gel(P,i) = RgX_div(V,W); e[i] = k; i++; }
1649 if (dW <= 0) break;
1650 V = W;
1651 }
1652 if (val) { gel(P,i) = pol_x[varn(f)]; e[i] = val; i++;}
1653 setlg(P,i);
1654 setlg(e,i); *ex = e; return P;
1655 }
1656
1657 GEN
fact_from_DDF(GEN fa,GEN e,long n)1658 fact_from_DDF(GEN fa, GEN e, long n)
1659 {
1660 GEN v,w, y = cgetg(3, t_MAT);
1661 long i,j,k, l = lg(fa);
1662
1663 v = cgetg(n+1,t_COL); gel(y,1) = v;
1664 w = cgetg(n+1,t_COL); gel(y,2) = w;
1665 for (k=i=1; i<l; i++)
1666 {
1667 GEN L = gel(fa,i), ex = utoipos(e[i]);
1668 long J = lg(L);
1669 for (j=1; j<J; j++,k++)
1670 {
1671 gel(v,k) = gcopy(gel(L,j));
1672 gel(w,k) = ex;
1673 }
1674 }
1675 return y;
1676 }
1677
1678 /* Factor x in Z[t]. Assume all factors have degree divisible by hint */
1679 GEN
factpol(GEN x,long hint)1680 factpol(GEN x, long hint)
1681 {
1682 pari_sp av = avma;
1683 GEN fa,ex,y;
1684 long n,i,l;
1685
1686 if (typ(x)!=t_POL) pari_err(notpoler,"factpol");
1687 if (!signe(x)) pari_err(zeropoler,"factpol");
1688
1689 fa = ZX_squff(Q_primpart(x), &ex);
1690 l = lg(fa); n = 0;
1691 for (i=1; i<l; i++)
1692 {
1693 gel(fa,i) = ZX_DDF(gel(fa,i), hint);
1694 n += lg(fa[i])-1;
1695 }
1696 y = fact_from_DDF(fa,ex,n);
1697 return gerepileupto(av, sort_factor(y, cmpii));
1698 }
1699
1700 GEN
nfrootsQ(GEN x)1701 nfrootsQ(GEN x)
1702 {
1703 pari_sp av = avma;
1704 GEN z, d;
1705 long val;
1706
1707 if (typ(x)!=t_POL) pari_err(notpoler,"nfrootsQ");
1708 if (!signe(x)) pari_err(zeropoler,"nfrootsQ");
1709 val = ZX_valuation(Q_primpart(x), &x);
1710 d = modulargcd(derivpol(x), x);
1711 if (degpol(d)) x = RgX_div(x, d);
1712 z = DDF(x, 1, 1);
1713 if (val) z = shallowconcat(z, gen_0);
1714 return gerepilecopy(av, z);
1715 }
1716
1717 /***********************************************************************/
1718 /** **/
1719 /** FACTORIZATION **/
1720 /** **/
1721 /***********************************************************************/
1722 #define LT 17
1723 #define assign_or_fail(x,y) {\
1724 if (y==NULL) y=x; else if (!gequal(x,y)) return 0;\
1725 }
1726 #define tsh 6
1727 #define typs(x,y) ((x << tsh) | y)
1728 #define typ1(x) (x >> tsh)
1729 #define typ2(x) (x & ((1<<tsh)-1))
1730
1731 static long
poltype(GEN x,GEN * ptp,GEN * ptpol,long * ptpa)1732 poltype(GEN x, GEN *ptp, GEN *ptpol, long *ptpa)
1733 {
1734 long t[LT]; /* code for 0,1,2,3,61,62,63,67,7,81,82,83,86,87,91,93,97 */
1735 long tx = typ(x),lx,i,j,s,pa=BIGINT;
1736 GEN pcx=NULL, p=NULL,pol=NULL,p1,p2;
1737
1738 if (is_scalar_t(tx))
1739 {
1740 if (tx == t_POLMOD) return 0;
1741 x = scalarpol(x,0);
1742 }
1743 for (i=2; i<LT; i++) t[i]=0; /* t[0..1] unused */
1744 lx = lg(x);
1745 for (i=2; i<lx; i++)
1746 {
1747 p1=gel(x,i);
1748 switch(typ(p1))
1749 {
1750 case t_INT: case t_FRAC:
1751 break;
1752 case t_REAL:
1753 s = precision(p1); if (s < pa) pa = s;
1754 t[2]=1; break;
1755 case t_INTMOD:
1756 assign_or_fail(gel(p1,1),p);
1757 t[3]=1; break;
1758 case t_COMPLEX:
1759 if (!pcx) pcx = mkpoln(3, gen_1,gen_0,gen_1); /* x^2 + 1 */
1760 for (j=1; j<=2; j++)
1761 {
1762 p2 = gel(p1,j);
1763 switch(typ(p2))
1764 {
1765 case t_INT: case t_FRAC:
1766 assign_or_fail(pcx,pol);
1767 t[4]=1; break;
1768 case t_REAL:
1769 s = precision(p2); if (s < pa) pa = s;
1770 t[5]=1; break;
1771 case t_INTMOD:
1772 assign_or_fail(gel(p2,1),p);
1773 assign_or_fail(pcx,pol);
1774 t[6]=1; break;
1775 case t_PADIC:
1776 s = precp(p2) + valp(p2); if (s < pa) pa = s;
1777 assign_or_fail(gel(p2,2),p);
1778 assign_or_fail(pcx,pol);
1779 t[7]=1; break;
1780 default: return 0;
1781 }
1782 }
1783 break;
1784 case t_PADIC:
1785 s = precp(p1) + valp(p1); if (s < pa) pa = s;
1786 assign_or_fail(gel(p1,2),p);
1787 t[8]=1; break;
1788 case t_QUAD:
1789 for (j=2; j<=3; j++)
1790 {
1791 p2 = gel(p1,j);
1792 switch(typ(p2))
1793 {
1794 case t_INT: case t_FRAC:
1795 assign_or_fail(gel(p1,1),pol);
1796 t[9]=1; break;
1797 case t_REAL:
1798 s = precision(p2); if (s < pa) pa = s;
1799 if (gsigne(discsr(gel(p1,1)))>0) t[10]=1; else t[12]=1;
1800 break;
1801 case t_INTMOD:
1802 assign_or_fail(gel(p2,1),p);
1803 assign_or_fail(gel(p1,1),pol);
1804 t[11]=1; break;
1805 case t_PADIC:
1806 s = precp(p2) + valp(p2); if (s < pa) pa = s;
1807 assign_or_fail(gel(p2,2),p);
1808 assign_or_fail(gel(p1,1),pol);
1809 t[13]=1; break;
1810 default: return 0;
1811 }
1812 }
1813 break;
1814 case t_POLMOD:
1815 assign_or_fail(gel(p1,1),pol);
1816 for (j=1; j<=2; j++)
1817 {
1818 GEN pbis = NULL, polbis = NULL;
1819 long pabis;
1820 switch(poltype(gel(p1,j),&pbis,&polbis,&pabis))
1821 {
1822 case t_INT: t[14]=1; break;
1823 case t_INTMOD: t[15]=1; break;
1824 case t_PADIC: t[16]=1; if (pabis<pa) pa=pabis; break;
1825 default: return 0;
1826 }
1827 if (pbis) assign_or_fail(pbis,p);
1828 if (polbis) assign_or_fail(polbis,pol);
1829 }
1830 break;
1831 default: return 0;
1832 }
1833 }
1834 if (t[5]||t[12])
1835 {
1836 if (t[3]||t[6]||t[7]||t[8]||t[11]||t[13]||t[14]||t[15]||t[16]) return 0;
1837 *ptpa=pa; return t_COMPLEX;
1838 }
1839 if (t[2]||t[10])
1840 {
1841 if (t[3]||t[6]||t[7]||t[8]||t[11]||t[13]||t[14]||t[15]||t[16]) return 0;
1842 *ptpa=pa; return t[4]?t_COMPLEX:t_REAL;
1843 }
1844 if (t[6]||t[11]||t[15])
1845 {
1846 *ptpol=pol; *ptp=p;
1847 i = t[15]? t_POLMOD: (t[11]? t_QUAD: t_COMPLEX);
1848 return typs(i, t_INTMOD);
1849 }
1850 if (t[7]||t[13]||t[16])
1851 {
1852 *ptpol=pol; *ptp=p; *ptpa=pa;
1853 i = t[16]? t_POLMOD: (t[13]? t_QUAD: t_COMPLEX);
1854 return typs(i, t_PADIC);
1855 }
1856 if (t[4]||t[9]||t[14])
1857 {
1858 *ptpol=pol;
1859 i = t[14]? t_POLMOD: (t[9]? t_QUAD: t_COMPLEX);
1860 return typs(i, t_INT);
1861 }
1862 if (t[3]) { *ptp=p; return t_INTMOD; }
1863 if (t[8]) { *ptp=p; *ptpa=pa; return t_PADIC; }
1864 return t_INT;
1865 }
1866 #undef LT
1867
1868 GEN
factor0(GEN x,long flag)1869 factor0(GEN x,long flag)
1870 {
1871 long tx=typ(x);
1872
1873 if (flag<0) return factor(x);
1874 if (is_matvec_t(tx)) return gboundfact(x,flag);
1875 if (tx==t_INT || tx==t_FRAC) return boundfact(x,flag);
1876 pari_err(talker,"partial factorization is not meaningful here");
1877 return NULL; /* not reached */
1878 }
1879
1880 GEN
concat_factor(GEN f,GEN g)1881 concat_factor(GEN f, GEN g)
1882 {
1883 if (lg(f) == 1) return g;
1884 if (lg(g) == 1) return f;
1885 return mkmat2(shallowconcat(gel(f,1), gel(g,1)),
1886 shallowconcat(gel(f,2), gel(g,2)));
1887 }
1888
1889 /* assume f and g coprime integer factorizations */
1890 GEN
merge_factor_i(GEN f,GEN g)1891 merge_factor_i(GEN f, GEN g)
1892 {
1893 if (lg(f) == 1) return g;
1894 if (lg(g) == 1) return f;
1895 return sort_factor_gen(concat_factor(f,g), cmpii);
1896 }
1897
1898 GEN
deg1_from_roots(GEN L,long v)1899 deg1_from_roots(GEN L, long v)
1900 {
1901 long i, l = lg(L);
1902 GEN z = cgetg(l,t_COL);
1903 for (i=1; i<l; i++)
1904 gel(z,i) = deg1pol_i(gen_1, gneg(gel(L,i)), v);
1905 return z;
1906 }
1907 GEN
roots_from_deg1(GEN x)1908 roots_from_deg1(GEN x)
1909 {
1910 long i,l = lg(x);
1911 GEN r = cgetg(l,t_VEC);
1912 for (i=1; i<l; i++) gel(r,i) = gneg(constant_term(gel(x,i)));
1913 return r;
1914 }
1915
1916 static GEN
gauss_factor_p(GEN p)1917 gauss_factor_p(GEN p)
1918 {
1919 GEN a, b; (void)cornacchia(gen_1, p, &a,&b);
1920 return mkcomplex(a, b);
1921 }
1922
1923 static GEN
gauss_primpart(GEN x,GEN * c)1924 gauss_primpart(GEN x, GEN *c)
1925 {
1926 GEN y, a = gel(x,1), b = gel(x,2), n = gcdii(a, b);
1927 *c = n; if (n == gen_1) return x;
1928 y = cgetg(3, t_COMPLEX);
1929 gel(y,1) = diviiexact(a, n);
1930 gel(y,2) = diviiexact(b, n); return y;
1931 }
1932
1933 static GEN
gauss_primpart_try(GEN x,GEN c)1934 gauss_primpart_try(GEN x, GEN c)
1935 {
1936 GEN r, y;
1937 if (typ(x) == t_INT)
1938 {
1939 y = dvmdii(x, c, &r); if (r != gen_0) return NULL;
1940 }
1941 else
1942 {
1943 GEN a = gel(x,1), b = gel(x,2); y = cgetg(3, t_COMPLEX);
1944 gel(y,1) = dvmdii(a, c, &r); if (r != gen_0) return NULL;
1945 gel(y,2) = dvmdii(b, c, &r); if (r != gen_0) return NULL;
1946 }
1947 return y;
1948 }
1949
1950 static int
gauss_cmp(GEN x,GEN y)1951 gauss_cmp(GEN x, GEN y)
1952 {
1953 int v;
1954 if (typ(x) != t_COMPLEX)
1955 return (typ(y) == t_COMPLEX)? -1: gcmp(x, y);
1956 if (typ(y) != t_COMPLEX) return 1;
1957 v = cmpii(gel(x,2), gel(y,2));
1958 if (v) return v;
1959 return gcmp(gel(x,1), gel(y,1));
1960 }
1961
1962 /* 0 or canonical representative in Z[i]^* / <i> (impose imag(x) >= 0) */
1963 static GEN
gauss_normal(GEN x)1964 gauss_normal(GEN x)
1965 {
1966 if (typ(x) != t_COMPLEX) return (signe(x) < 0)? absi(x): x;
1967 if (signe(x[1]) < 0) x = gneg(x);
1968 if (signe(x[2]) < 0) x = mulcxI(x);
1969 return x;
1970 }
1971
1972 static GEN
Ipow(long e)1973 Ipow(long e) {
1974 switch(e & 3)
1975 {
1976 case 1: return gi;
1977 case 2: return gen_m1;
1978 case 3: return pureimag(gen_m1);
1979 }
1980 return gen_1;
1981 }
1982
1983 static GEN
gauss_factor(GEN x)1984 gauss_factor(GEN x)
1985 {
1986 pari_sp av = avma;
1987 GEN a = gel(x,1), b = gel(x,2), d = gen_1, n, y, fa, P, E, P2, E2;
1988 long t1 = typ(a);
1989 long t2 = typ(b), i, j, l, exp = 0;
1990 if (t1 == t_FRAC) d = gel(a,2);
1991 if (t2 == t_FRAC) d = lcmii(d, gel(b,2));
1992 if (d == gen_1) y = x;
1993 else
1994 {
1995 y = gmul(x, d);
1996 a = gel(y,1); t1 = typ(a);
1997 b = gel(y,2); t2 = typ(b);
1998 }
1999 if (t1 != t_INT || t2 != t_INT) return NULL;
2000 y = gauss_primpart(y, &n);
2001 fa = factor(cxnorm(y));
2002 P = gel(fa,1);
2003 E = gel(fa,2); l = lg(P);
2004 P2 = cgetg(l, t_COL);
2005 E2 = cgetg(l, t_COL);
2006 for (j = 1, i = l-1; i > 0; i--) /* remove largest factors first */
2007 {
2008 GEN p = gel(P,i), w, w2, t, we, pe;
2009 long v, e = itos(gel(E,i));
2010 int is2 = equaliu(p, 2);
2011 w = is2? mkcomplex(gen_1,gen_1): gauss_factor_p(p);
2012 w2 = gauss_normal( gconj(w) );
2013 /* w * w2 * I^3 = p, w2 = gconj(w) * I */
2014 pe = powiu(p, e);
2015 we = gpowgs(w, e);
2016 t = gauss_primpart_try( gmul(y, gconj(we)), pe );
2017 if (t) y = t; /* y /= w^e */
2018 else {
2019 /* y /= conj(w)^e, should be y /= w2^e */
2020 y = gauss_primpart_try( gmul(y, we), pe );
2021 swap(w, w2); exp += 3 * e;
2022 }
2023 gel(P,i) = w;
2024 v = Z_pvalrem(n, p, &n);
2025 if (v) {
2026 exp += 3*v;
2027 if (is2) v <<= 1; /* 2 = w^2 I^3 */
2028 else {
2029 gel(P2,j) = w2;
2030 gel(E2,j) = utoipos(v); j++;
2031 }
2032 gel(E,i) = stoi(e + v);
2033 }
2034 v = Z_pvalrem(d, p, &d);
2035 if (v) {
2036 exp -= 3*v;
2037 if (is2) v <<= 1; /* 2 is ramified */
2038 else {
2039 gel(P2,j) = w2;
2040 gel(E2,j) = utoineg(v); j++;
2041 }
2042 gel(E,i) = stoi(e - v);
2043 }
2044 exp &= 3;
2045 }
2046 if (j > 1) {
2047 setlg(P2, j);
2048 setlg(E2, j);
2049 fa = concat_factor(fa, mkmat2(P2,E2));
2050 }
2051 if (!is_pm1(n) || !is_pm1(d))
2052 {
2053 GEN Fa = factor(gdiv(n, d));
2054 P = gel(Fa,1); l = lg(P);
2055 E = gel(Fa,2);
2056 for (i = j = 1; i < l; i++)
2057 {
2058 GEN w, p = gel(P,i);
2059 long e;
2060 int is2;
2061 switch(mod4(p))
2062 {
2063 case 3: continue;
2064 case 2: is2 = 1; break;
2065 default:is2 = 0; break;
2066 }
2067 e = itos(gel(E,i));
2068 w = is2? mkcomplex(gen_1,gen_1): gauss_factor_p(p);
2069 gel(P,i) = w;
2070 if (is2)
2071 gel(E,i) = stoi(e << 1);
2072 else
2073 {
2074 P = shallowconcat(P, gauss_normal( gconj(w) ));
2075 E = shallowconcat(E, gel(E,i));
2076 }
2077 exp += 3*e;
2078 exp &= 3;
2079 }
2080 gel(Fa,1) = P;
2081 gel(Fa,2) = E;
2082 fa = concat_factor(fa, Fa);
2083 }
2084 fa = sort_factor_gen(fa, &gauss_cmp);
2085
2086 y = gmul(y, Ipow(exp));
2087 if (!gcmp1(y)) {
2088 gel(fa,1) = shallowconcat(mkcol(y), gel(fa,1));
2089 gel(fa,2) = shallowconcat(gen_1, gel(fa,2));
2090 }
2091 return gerepilecopy(av, fa);
2092 }
2093
2094 GEN
factor(GEN x)2095 factor(GEN x)
2096 {
2097 long tx=typ(x), lx, i, j, pa, v, r1;
2098 pari_sp av, tetpil;
2099 GEN y,p,p1,p2,p3,p4,p5,pol;
2100
2101 if (is_matvec_t(tx))
2102 {
2103 lx=lg(x); y=cgetg(lx,tx);
2104 for (i=1; i<lx; i++) gel(y,i) = factor(gel(x,i));
2105 return y;
2106 }
2107 if (gcmp0(x))
2108 {
2109 y = cgetg(3,t_MAT);
2110 gel(y,1) = mkcolcopy(x);
2111 gel(y,2) = mkcol(gen_1); return y;
2112 }
2113 av = avma;
2114 switch(tx)
2115 {
2116 case t_INT: return Z_factor(x);
2117
2118 case t_FRAC:
2119 p1 = Z_factor(gel(x,1));
2120 p2 = Z_factor(gel(x,2)); gel(p2,2) = gneg_i(gel(p2,2));
2121 return gerepilecopy(av, merge_factor_i(p1,p2));
2122
2123 case t_POL:
2124 tx=poltype(x,&p,&pol,&pa);
2125 switch(tx)
2126 {
2127 case 0: pari_err(impl,"factor for general polynomials");
2128 case t_INT: return factpol(x,1);
2129 case t_INTMOD: return factmod(x,p);
2130
2131 case t_COMPLEX: y=cgetg(3,t_MAT); lx=lg(x)-2; v=varn(x);
2132 av = avma; p1 = roots(x,pa); tetpil = avma;
2133 p2 = deg1_from_roots(p1, v);
2134 gel(y,1) = gerepile(av,tetpil,p2);
2135 p3=cgetg(lx,t_COL); for (i=1; i<lx; i++) gel(p3,i) = gen_1;
2136 gel(y,2) = p3; return y;
2137
2138 case t_REAL: y=cgetg(3,t_MAT); lx=lg(x)-2; v=varn(x);
2139 av=avma; p1=cleanroots(x,pa); tetpil=avma;
2140 for(r1=1; r1<lx; r1++)
2141 if (typ(p1[r1]) == t_COMPLEX) break;
2142 lx=(r1+lx)>>1; p2=cgetg(lx,t_COL);
2143 for(i=1; i<r1; i++)
2144 gel(p2,i) = deg1pol_i(gen_1, negr(gel(p1,i)), v);
2145 for( ; i<lx; i++)
2146 {
2147 GEN a = gel(p1,2*i-r1);
2148 p = cgetg(5, t_POL); gel(p2,i) = p;
2149 p[1] = x[1];
2150 gel(p,2) = gnorm(a);
2151 gel(p,3) = gmul2n(gel(a,1),1); setsigne(p[3],-signe(p[3]));
2152 gel(p,4) = gen_1;
2153 }
2154 gel(y,1) = gerepile(av,tetpil,p2);
2155 p3=cgetg(lx,t_COL); for (i=1; i<lx; i++) gel(p3,i) = gen_1;
2156 gel(y,2) = p3; return y;
2157
2158 case t_PADIC: return factorpadic4(x,p,pa);
2159
2160 default:
2161 {
2162 long killv;
2163 x = shallowcopy(x); lx=lg(x);
2164 pol = shallowcopy(pol);
2165 v = manage_var(manage_var_max_avail,NULL);
2166 for(i=2; i<lx; i++)
2167 {
2168 p1=gel(x,i);
2169 switch(typ(p1))
2170 {
2171 case t_QUAD: p1++;
2172 case t_COMPLEX:
2173 p2 = cgetg(3, t_POLMOD); gel(x,i) = p2;
2174 gel(p2,1) = pol;
2175 gel(p2,2) = deg1pol_i(gel(p1,2), gel(p1,1), v);
2176 }
2177 }
2178 killv = (avma != (pari_sp)pol);
2179 if (killv) setvarn(pol, fetch_var());
2180 switch (typ2(tx))
2181 {
2182 case t_INT: p1 = polfnf(x,pol); break;
2183 case t_INTMOD: p1 = factorff(x,p,pol); break;
2184 default: pari_err(impl,"factor of general polynomial");
2185 return NULL; /* not reached */
2186 }
2187 switch (typ1(tx))
2188 {
2189 case t_POLMOD:
2190 if (killv) (void)delete_var();
2191 return gerepileupto(av,p1);
2192 case t_COMPLEX: p5 = gi; break;
2193 case t_QUAD: p5=cgetg(4,t_QUAD);
2194 gel(p5,1) = pol; gel(p5,2) = gen_0; gel(p5,3) = gen_1;
2195 break;
2196 default: pari_err(impl,"factor of general polynomial");
2197 return NULL; /* not reached */
2198 }
2199 p2=gel(p1,1);
2200 for(i=1; i<lg(p2); i++)
2201 {
2202 p3=gel(p2,i);
2203 for(j=2; j<lg(p3); j++)
2204 {
2205 p4=gel(p3,j);
2206 if(typ(p4)==t_POLMOD) gel(p3,j) = gsubst(gel(p4,2),v,p5);
2207 }
2208 }
2209 if (killv) (void)delete_var();
2210 return gerepilecopy(av, p1);
2211 }
2212 }
2213 case t_RFRAC:
2214 p1 = factor(gel(x,1));
2215 p2 = factor(gel(x,2)); gel(p2,2) = gneg_i(gel(p2,2));
2216 return gerepilecopy(av, concat_factor(p1,p2));
2217
2218 case t_COMPLEX:
2219 y = gauss_factor(x);
2220 if (y) return y;
2221 }
2222 pari_err(talker,"can't factor %Z",x);
2223 return NULL; /* not reached */
2224 }
2225 #undef typ1
2226 #undef typ2
2227 #undef typs
2228 #undef tsh
2229
2230 /* assume n > 0. Compute x^n using left-right binary powering */
2231 GEN
leftright_pow_u_fold(GEN x,ulong n,void * data,GEN (* sqr)(void *,GEN),GEN (* msqr)(void *,GEN))2232 leftright_pow_u_fold(GEN x, ulong n, void *data, GEN (*sqr)(void*,GEN),
2233 GEN (*msqr)(void*,GEN))
2234 {
2235 GEN y;
2236 long m, j;
2237
2238 if (n == 1) return gcopy(x);
2239
2240 m = (long)n; j = 1+bfffo(m);
2241 y = x;
2242
2243 /* normalize, i.e set highest bit to 1 (we know m != 0) */
2244 m<<=j; j = BITS_IN_LONG-j;
2245 /* first bit is now implicit */
2246 for (; j; m<<=1,j--)
2247 {
2248 if (m < 0) y = msqr(data,y); /* first bit set: multiply by base */
2249 else y = sqr(data,y);
2250 }
2251 return y;
2252 }
2253
2254
2255 /* assume n != 0, t_INT. Compute x^|n| using left-right binary powering */
2256 GEN
leftright_pow_fold(GEN x,GEN n,void * data,GEN (* sqr)(void *,GEN),GEN (* msqr)(void *,GEN))2257 leftright_pow_fold(GEN x, GEN n, void *data, GEN (*sqr)(void*,GEN),
2258 GEN (*msqr)(void*,GEN))
2259 {
2260 long ln = lgefint(n);
2261 if (ln == 3) return leftright_pow_u_fold(x, n[2], data, sqr, msqr);
2262 else
2263 {
2264 GEN nd = int_MSW(n), y = x;
2265 long i, m = *nd, j = 1+bfffo((ulong)m);
2266 pari_sp av = avma, lim = stack_lim(av, 1);
2267
2268 /* normalize, i.e set highest bit to 1 (we know m != 0) */
2269 m<<=j; j = BITS_IN_LONG-j;
2270 /* first bit is now implicit */
2271 for (i=ln-2;;)
2272 {
2273 for (; j; m<<=1,j--)
2274 {
2275 if (m < 0) y = msqr(data,y); /* first bit set: multiply by base */
2276 else y = sqr(data,y);
2277 if (low_stack(lim, stack_lim(av,1)))
2278 {
2279 if (DEBUGMEM>1) pari_warn(warnmem,"leftright_pow");
2280 y = gerepilecopy(av, y);
2281 }
2282 }
2283 if (--i == 0) return y;
2284 nd=int_precW(nd);
2285 m = *nd; j = BITS_IN_LONG;
2286 }
2287 }
2288 }
2289
2290 struct leftright_fold
2291 {
2292 void *data;
2293 GEN x;
2294 GEN (*sqr)(void*,GEN);
2295 GEN (*mul)(void*,GEN,GEN);
2296 };
2297
2298 static GEN
leftright_sqr(void * data,GEN y)2299 leftright_sqr(void* data, GEN y)
2300 {
2301 struct leftright_fold *d=(struct leftright_fold*) data;
2302 return d->sqr(d->data,y);
2303 }
2304
2305 static GEN
leftright_msqr(void * data,GEN y)2306 leftright_msqr(void* data, GEN y)
2307 {
2308 struct leftright_fold *d=(struct leftright_fold*) data;
2309 return d->mul(d->data,d->sqr(d->data,y),d->x);
2310 }
2311
2312 GEN
leftright_pow(GEN x,GEN n,void * data,GEN (* sqr)(void *,GEN),GEN (* mul)(void *,GEN,GEN))2313 leftright_pow(GEN x, GEN n, void *data, GEN (*sqr)(void*,GEN),
2314 GEN (*mul)(void*,GEN,GEN))
2315 {
2316 struct leftright_fold d;
2317 d.data=data; d.mul=mul; d.sqr=sqr; d.x=x;
2318 return leftright_pow_fold(x, n, (void *)&d, leftright_sqr, leftright_msqr);
2319 }
2320
2321 GEN
leftright_pow_u(GEN x,ulong n,void * data,GEN (* sqr)(void *,GEN),GEN (* mul)(void *,GEN,GEN))2322 leftright_pow_u(GEN x, ulong n, void *data, GEN (*sqr)(void*,GEN),
2323 GEN (*mul)(void*,GEN,GEN))
2324 {
2325 struct leftright_fold d;
2326 d.data=data; d.mul=mul; d.sqr=sqr; d.x=x;
2327 return leftright_pow_u_fold(x, n, (void *)&d, leftright_sqr, leftright_msqr);
2328 }
2329
2330 GEN
divide_conquer_assoc(GEN x,GEN (* mul)(void *,GEN,GEN),void * data)2331 divide_conquer_assoc(GEN x, GEN (*mul)(void *,GEN,GEN),void *data)
2332 {
2333 pari_sp ltop, lim;
2334 long i,k,lx = lg(x);
2335
2336 if (lx == 1) return gen_1;
2337 if (lx == 2) return gcopy(gel(x,1));
2338 x = shallowcopy(x); k = lx;
2339 ltop=avma; lim = stack_lim(ltop,1);
2340 while (k > 2)
2341 {
2342 if (DEBUGLEVEL>7)
2343 fprintferr("prod: remaining objects %ld\n",k-1);
2344 lx = k; k = 1;
2345 for (i=1; i<lx-1; i+=2)
2346 gel(x,k++) = mul(data,gel(x,i),gel(x,i+1));
2347 if (i < lx) x[k++] = x[i];
2348 if (low_stack(lim,stack_lim(av,1)))
2349 gerepilecoeffs(ltop,x+1,k-1);
2350 }
2351 return gel(x,1);
2352 }
2353
2354 static GEN
_domul(void * data,GEN x,GEN y)2355 _domul(void *data, GEN x, GEN y)
2356 {
2357 GEN (*mul)(GEN,GEN)=(GEN (*)(GEN,GEN)) data;
2358 return mul(x,y);
2359 }
2360
2361 GEN
divide_conquer_prod(GEN x,GEN (* mul)(GEN,GEN))2362 divide_conquer_prod(GEN x, GEN (*mul)(GEN,GEN))
2363 {
2364 return divide_conquer_assoc(x, _domul, (void *)mul);
2365 }
2366
2367 static GEN
idmulred(void * nf,GEN x,GEN y)2368 idmulred(void *nf, GEN x, GEN y) { return idealmulred((GEN) nf, x, y, 0); }
2369 static GEN
idpowred(void * nf,GEN x,GEN n)2370 idpowred(void *nf, GEN x, GEN n) { return idealpowred((GEN) nf, x, n, 0); }
2371 static GEN
idmul(void * nf,GEN x,GEN y)2372 idmul(void *nf, GEN x, GEN y) { return idealmul((GEN) nf, x, y); }
2373 static GEN
idpow(void * nf,GEN x,GEN n)2374 idpow(void *nf, GEN x, GEN n) { return idealpow((GEN) nf, x, n); }
2375 static GEN
eltmul(void * nf,GEN x,GEN y)2376 eltmul(void *nf, GEN x, GEN y) { return element_mul((GEN) nf, x, y); }
2377 static GEN
eltpow(void * nf,GEN x,GEN n)2378 eltpow(void *nf, GEN x, GEN n) { return element_pow((GEN) nf, x, n); }
2379
2380 #if 0
2381 static GEN
2382 ellmul(void *ell, GEN x, GEN y) { return addell((GEN) ell, x, y); }
2383 static GEN
2384 ellpow(void *GEN x, GEN n) { return powell((GEN) ell, x, n); }
2385 #endif
2386
2387 GEN
factorback_aux(GEN fa,GEN e,GEN (* _mul)(void *,GEN,GEN),GEN (* _pow)(void *,GEN,GEN),void * data)2388 factorback_aux(GEN fa, GEN e, GEN (*_mul)(void*,GEN,GEN), GEN (*_pow)(void*,GEN,GEN), void *data)
2389 {
2390 pari_sp av = avma;
2391 long k,l,lx,t = typ(fa);
2392 GEN p,x;
2393
2394 if (e) /* supplied vector of exponents */
2395 p = fa;
2396 else /* genuine factorization */
2397 {
2398 if (t == t_MAT) {
2399 l = lg(fa);
2400 if (l == 1) return gen_1;
2401 if (l != 3) pari_err(talker,"not a factorisation in factorback");
2402 } else {
2403 if (!is_vec_t(t)) pari_err(talker,"not a factorisation in factorback");
2404 return gerepileupto(av, divide_conquer_assoc(fa, _mul,data));
2405 }
2406 p = gel(fa,1);
2407 e = gel(fa,2);
2408 }
2409 lx = lg(p);
2410 t = t_INT; /* dummy */
2411 /* check whether e is an integral vector of correct length */
2412 if (is_vec_t(typ(e)) && lx == lg(e))
2413 {
2414 for (k=1; k<lx; k++)
2415 if (typ(e[k]) != t_INT) break;
2416 if (k == lx) t = t_MAT;
2417 }
2418 if (t != t_MAT) pari_err(talker,"not a factorisation in factorback");
2419 if (lx == 1) return gen_1;
2420 x = cgetg(lx,t_VEC);
2421 for (l=1,k=1; k<lx; k++)
2422 if (signe(e[k]))
2423 gel(x,l++) = _pow(data,gel(p,k),gel(e,k));
2424 setlg(x,l);
2425 return gerepileupto(av, divide_conquer_assoc(x, _mul,data));
2426 }
2427
_agmul(void * a,GEN x,GEN y)2428 static GEN _agmul(void *a, GEN x, GEN y) { return gmul(x,y);}
_apowgi(void * a,GEN x,GEN y)2429 static GEN _apowgi(void *a, GEN x, GEN y) { return powgi(x,y);}
2430
2431 GEN
factorback_i(GEN fa,GEN e,GEN OBJ,int red)2432 factorback_i(GEN fa, GEN e, GEN OBJ, int red)
2433 {
2434 if (!OBJ)
2435 {
2436 if (e) {
2437 OBJ = checknf_i(e); if (OBJ) e = NULL;
2438 }
2439 if (!OBJ) return factorback_aux(fa, e, &_agmul, &_apowgi, NULL);
2440 }
2441 if (red) return factorback_aux(fa, e, &idmulred, &idpowred, OBJ);
2442 else return factorback_aux(fa, e, &idmul, &idpow, OBJ);
2443 }
2444
2445 GEN
factorbackelt(GEN fa,GEN e,GEN nf)2446 factorbackelt(GEN fa, GEN e, GEN nf)
2447 {
2448 if (!nf && e && lg(e) > 1 && typ(e[1]) != t_INT) { nf = e; e = NULL; }
2449 if (!nf) pari_err(talker, "missing nf in factorbackelt");
2450
2451 nf = checknf(nf);
2452 return factorback_aux(fa, e, &eltmul, &eltpow, nf);
2453 }
2454
2455 GEN
factorback0(GEN fa,GEN e,GEN nf)2456 factorback0(GEN fa, GEN e, GEN nf)
2457 {
2458 return factorback_i(fa,e,nf,0);
2459 }
2460
2461 GEN
factorback(GEN fa,GEN nf)2462 factorback(GEN fa, GEN nf)
2463 {
2464 return factorback_i(fa,nf,NULL,0);
2465 }
2466
2467 GEN
gisirreducible(GEN x)2468 gisirreducible(GEN x)
2469 {
2470 long tx = typ(x), l, i;
2471 pari_sp av=avma;
2472 GEN y;
2473
2474 if (is_matvec_t(tx))
2475 {
2476 l=lg(x); y=cgetg(l,tx);
2477 for (i=1; i<l; i++) gel(y,i) = gisirreducible(gel(x,i));
2478 return y;
2479 }
2480 if (is_intreal_t(tx) || tx == t_FRAC) return gen_0;
2481 if (tx!=t_POL) pari_err(notpoler,"gisirreducible");
2482 l=lg(x); if (l<=3) return gen_0;
2483 y=factor(x); avma=av;
2484 return (lg(gcoeff(y,1,1))==l)?gen_1:gen_0;
2485 }
2486
2487 /*******************************************************************/
2488 /* */
2489 /* GENERIC GCD */
2490 /* */
2491 /*******************************************************************/
2492 GEN
gcd0(GEN x,GEN y,long flag)2493 gcd0(GEN x, GEN y, long flag)
2494 {
2495 if (!y) return content(x);
2496 switch(flag)
2497 {
2498 case 0: return ggcd(x,y);
2499 case 1: return modulargcd(x,y);
2500 case 2: return srgcd(x,y);
2501 default: pari_err(flagerr,"gcd");
2502 }
2503 return NULL; /* not reached */
2504 }
2505
2506 /* x is a COMPLEX or a QUAD */
2507 static GEN
triv_cont_gcd(GEN x,GEN y)2508 triv_cont_gcd(GEN x, GEN y)
2509 {
2510 pari_sp av = avma, tetpil;
2511 GEN p1 = (typ(x)==t_COMPLEX)? ggcd(gel(x,1),gel(x,2))
2512 : ggcd(gel(x,2),gel(x,3));
2513 tetpil=avma; return gerepile(av,tetpil,ggcd(p1,y));
2514 }
2515
2516 /* y is a PADIC, x a rational number or an INTMOD */
2517 static GEN
padic_gcd(GEN x,GEN y)2518 padic_gcd(GEN x, GEN y)
2519 {
2520 GEN p = gel(y,2);
2521 long v = ggval(x,p), w = valp(y);
2522 if (w < v) v = w;
2523 return gpowgs(p, v);
2524 }
2525
2526 /* x,y in Z[i], at least one of which is t_COMPLEX */
2527 static GEN
gauss_gcd(GEN x,GEN y)2528 gauss_gcd(GEN x, GEN y)
2529 {
2530 pari_sp av = avma;
2531 GEN dx, dy;
2532 dx = denom(x); x = gmul(x, dx);
2533 dy = denom(y); y = gmul(y, dy);
2534 while (!gcmp0(y))
2535 {
2536 GEN z = gsub(x, gmul(ground(gdiv(x,y)), y));
2537 x = y; y = z;
2538 }
2539 x = gauss_normal(x);
2540 if (typ(x) == t_COMPLEX)
2541 {
2542 if (gcmp0(gel(x,2))) x = gel(x,1);
2543 else if (gcmp0(gel(x,1))) x = gel(x,2);
2544 }
2545 return gerepileupto(av, gdiv(x, lcmii(dx, dy)));
2546 }
2547
2548 static int
c_is_rational(GEN x)2549 c_is_rational(GEN x)
2550 {
2551 return (is_rational(gel(x,1)) && is_rational(gel(x,2)));
2552 }
2553 static GEN
c_zero_gcd(GEN c)2554 c_zero_gcd(GEN c)
2555 {
2556 GEN x = gel(c,1), y = gel(c,2);
2557 long tx = typ(x), ty = typ(y);
2558 if (tx == t_REAL || ty == t_REAL) return gen_1;
2559 if (tx == t_PADIC || tx == t_INTMOD
2560 || ty == t_PADIC || ty == t_INTMOD) return ggcd(x, y);
2561 return gauss_gcd(c, gen_0);
2562 }
2563
2564 /* y == 0 */
2565 static GEN
zero_gcd(GEN x,long tx)2566 zero_gcd(GEN x, long tx)
2567 {
2568 pari_sp av;
2569 switch(tx)
2570 {
2571 case t_INT: return absi(x);
2572 case t_FRAC: return gabs(x,0);
2573 case t_COMPLEX: return c_zero_gcd(x);
2574 case t_REAL: return gen_1;
2575 case t_PADIC: return gpowgs(gel(x,2), valp(x));
2576 case t_SER: return monomial(gen_1, valp(x), varn(x));
2577 case t_POLMOD: {
2578 GEN d = gel(x,2);
2579 if (typ(d) == t_POL && varn(d) == varn(gel(x,1))) return content(d);
2580 return isinexact(d)? zero_gcd(d, typ(d)): gcopy(d);
2581 }
2582 case t_POL:
2583 if (!isinexact(x)) break;
2584 av = avma;
2585 return gerepileupto(av,
2586 monomialcopy(content(x), polvaluation(x,NULL), varn(x))
2587 );
2588
2589 case t_RFRAC:
2590 if (!isinexact(x)) break;
2591 av = avma;
2592 return gerepileupto(av,
2593 gdiv(zero_gcd(gel(x,1), typ(gel(x,1))), gel(x,2))
2594 );
2595 }
2596 return gcopy(x);
2597 }
2598
2599 /* tx = t_RFRAC, y considered as constant */
2600 static GEN
cont_gcd_rfrac(GEN x,GEN y)2601 cont_gcd_rfrac(GEN x, GEN y)
2602 {
2603 pari_sp av = avma;
2604 GEN cx; x = primitive_part(x, &cx);
2605 return gerepileupto(av, gred_rfrac_simple(ggcd(cx? cx: gen_1, y), gel(x,2)));
2606 }
2607 /* !is_const(tx), tx != t_RFRAC, y considered as constant */
2608 static GEN
cont_gcd_gen(GEN x,GEN y)2609 cont_gcd_gen(GEN x, GEN y)
2610 {
2611 pari_sp av = avma;
2612 return gerepileupto(av, ggcd(content(x),y));
2613 }
2614 /* !is_const(tx), y considered as constant */
2615 static GEN
cont_gcd(GEN x,long tx,GEN y)2616 cont_gcd(GEN x, long tx, GEN y)
2617 {
2618 return (tx == t_RFRAC)? cont_gcd_rfrac(x, y): cont_gcd_gen(x, y);
2619 }
2620
2621 GEN
ggcd(GEN x,GEN y)2622 ggcd(GEN x, GEN y)
2623 {
2624 long l, i, vx, vy, tx = typ(x), ty = typ(y);
2625 pari_sp av, tetpil;
2626 GEN p1,z;
2627
2628 if (is_noncalc_t(tx) || is_noncalc_t(ty)) pari_err(operf,"g",x,y);
2629 if (is_matvec_t(ty))
2630 {
2631 l = lg(y); z = cgetg(l,ty);
2632 for (i=1; i<l; i++) gel(z,i) = ggcd(x,gel(y,i));
2633 return z;
2634 }
2635 if (is_matvec_t(tx))
2636 {
2637 l = lg(x); z = cgetg(l,tx);
2638 for (i=1; i<l; i++) gel(z,i) = ggcd(gel(x,i),y);
2639 return z;
2640 }
2641 if (tx>ty) { swap(x,y); lswap(tx,ty); }
2642 /* tx <= ty */
2643 if (isexactzero(x)) return zero_gcd(y, ty);
2644 if (isexactzero(y)) return zero_gcd(x, tx);
2645 if (is_const_t(tx))
2646 {
2647 if (ty == tx) switch(tx)
2648 {
2649 case t_INT:
2650 return gcdii(x,y);
2651
2652 case t_INTMOD: z=cgetg(3,t_INTMOD);
2653 if (equalii(gel(x,1),gel(y,1)))
2654 gel(z,1) = gcopy(gel(x,1));
2655 else
2656 gel(z,1) = gcdii(gel(x,1),gel(y,1));
2657 if (gcmp1(gel(z,1))) gel(z,2) = gen_0;
2658 else
2659 {
2660 av = avma; p1 = gcdii(gel(z,1),gel(x,2));
2661 if (!is_pm1(p1)) p1 = gerepileuptoint(av, gcdii(p1,gel(y,2)));
2662 gel(z,2) = p1;
2663 }
2664 return z;
2665
2666 case t_FRAC: z=cgetg(3,t_FRAC);
2667 gel(z,1) = gcdii(gel(x,1), gel(y,1));
2668 gel(z,2) = lcmii(gel(x,2), gel(y,2));
2669 return z;
2670
2671 case t_COMPLEX:
2672 if (c_is_rational(x) && c_is_rational(y)) return gauss_gcd(x,y);
2673 return triv_cont_gcd(y,x);
2674
2675 case t_PADIC:
2676 if (!equalii(gel(x,2),gel(y,2))) return gen_1;
2677 vx = valp(x);
2678 vy = valp(y); return gpowgs(gel(y,2), min(vy,vx));
2679
2680 case t_QUAD:
2681 av=avma; p1=gdiv(x,y);
2682 if (gcmp0(gel(p1,3)))
2683 {
2684 p1=gel(p1,2);
2685 if (typ(p1)==t_INT) { avma=av; return gcopy(y); }
2686 tetpil=avma; return gerepile(av,tetpil, gdiv(y,gel(p1,2)));
2687 }
2688 if (typ(p1[2])==t_INT && typ(p1[3])==t_INT) {avma=av; return gcopy(y);}
2689 p1 = ginv(p1); avma=av;
2690 if (typ(p1[2])==t_INT && typ(p1[3])==t_INT) return gcopy(x);
2691 return triv_cont_gcd(y,x);
2692
2693 default: return gen_1; /* t_REAL */
2694 }
2695 if (is_const_t(ty)) switch(tx)
2696 {
2697 case t_INT:
2698 switch(ty)
2699 {
2700 case t_INTMOD: z = cgetg(3,t_INTMOD);
2701 gel(z,1) = gcopy(gel(y,1)); av = avma;
2702 p1 = gcdii(gel(y,1),gel(y,2));
2703 if (!is_pm1(p1)) p1 = gerepileuptoint(av, gcdii(x,p1));
2704 gel(z,2) = p1; return z;
2705
2706 case t_FRAC: z = cgetg(3,t_FRAC);
2707 gel(z,1) = gcdii(x,gel(y,1));
2708 gel(z,2) = icopy(gel(y,2)); return z;
2709
2710 case t_COMPLEX:
2711 if (c_is_rational(y)) return gauss_gcd(x,y);
2712 case t_QUAD:
2713 return triv_cont_gcd(y,x);
2714
2715 case t_PADIC:
2716 return padic_gcd(x,y);
2717
2718 default: return gen_1; /* t_REAL */
2719 }
2720
2721 case t_INTMOD:
2722 switch(ty)
2723 {
2724 case t_FRAC:
2725 av = avma; p1=gcdii(gel(x,1),gel(y,2)); avma = av;
2726 if (!is_pm1(p1)) pari_err(operi,"g",x,y);
2727 return ggcd(gel(y,1), x);
2728
2729 case t_COMPLEX: case t_QUAD:
2730 return triv_cont_gcd(y,x);
2731
2732 case t_PADIC:
2733 return padic_gcd(x,y);
2734 }
2735
2736 case t_FRAC:
2737 switch(ty)
2738 {
2739 case t_COMPLEX:
2740 if (c_is_rational(y)) return gauss_gcd(x,y);
2741 case t_QUAD:
2742 return triv_cont_gcd(y,x);
2743
2744 case t_PADIC:
2745 return padic_gcd(x,y);
2746 }
2747
2748 case t_COMPLEX: /* ty = PADIC or QUAD */
2749 return triv_cont_gcd(x,y);
2750
2751 case t_PADIC: /* ty = QUAD */
2752 return triv_cont_gcd(y,x);
2753
2754 default: return gen_1; /* tx = t_REAL */
2755 }
2756 return cont_gcd(y,ty, x);
2757 }
2758
2759 if (tx == t_POLMOD)
2760 {
2761 if (ty == t_POLMOD)
2762 {
2763 z = cgetg(3,t_POLMOD);
2764 if (gequal(gel(x,1),gel(y,1)))
2765 gel(z,1) = gcopy(gel(x,1));
2766 else
2767 gel(z,1) = ggcd(gel(x,1),gel(y,1));
2768 if (degpol(z[1])<=0) gel(z,2) = gen_0;
2769 else
2770 {
2771 GEN X, Y, d;
2772 av = avma; X = gel(x,2); Y = gel(y,2);
2773 d = ggcd(content(X), content(Y));
2774 if (!gcmp1(d)) { X = gdiv(X,d); Y = gdiv(Y,d); }
2775 p1 = ggcd(gel(z,1), X);
2776 gel(z,2) = gerepileupto(av, gmul(d, ggcd(p1, Y)));
2777 }
2778 return z;
2779 }
2780 vx = varn(x[1]);
2781 switch(ty)
2782 {
2783 case t_POL:
2784 vy = varn(y);
2785 if (varncmp(vy,vx) < 0) return cont_gcd_gen(y, x);
2786 z = cgetg(3,t_POLMOD);
2787 gel(z,1) = gcopy(gel(x,1));
2788 av = avma; p1 = ggcd(gel(x,1),gel(x,2));
2789 gel(z,2) = gerepileupto(av, ggcd(p1,y));
2790 return z;
2791
2792 case t_RFRAC:
2793 vy = varn(y[2]);
2794 if (varncmp(vy,vx) < 0) return cont_gcd_rfrac(y, x);
2795 av = avma;
2796 p1 = ggcd(gel(x,1),gel(y,2));
2797 if (degpol(p1)) pari_err(operi,"g",x,y);
2798 avma = av; return gdiv(ggcd(gel(y,1),x), content(gel(y,2)));
2799 }
2800 }
2801
2802 vx = gvar(x);
2803 vy = gvar(y);
2804 if (varncmp(vy, vx) < 0) return cont_gcd(y,ty, x);
2805 if (varncmp(vy, vx) > 0) return cont_gcd(x,tx, y);
2806
2807 /* same main variable */
2808 switch(tx)
2809 {
2810 case t_POL:
2811 switch(ty)
2812 {
2813 case t_POL: return srgcd(x,y);
2814 case t_SER:
2815 z = ggcd(content(x), content(y));
2816 return monomialcopy(z, min(valp(y),gval(x,vx)), vx);
2817 case t_RFRAC: return cont_gcd_rfrac(y, x);
2818 }
2819 break;
2820
2821 case t_SER:
2822 z = ggcd(content(x), content(y));
2823 switch(ty)
2824 {
2825 case t_SER: return monomialcopy(z, min(valp(x),valp(y)), vx);
2826 case t_RFRAC: return monomialcopy(z, min(valp(x),gval(y,vx)), vx);
2827 }
2828 break;
2829
2830 case t_RFRAC: z=cgetg(3,t_RFRAC);
2831 if (ty != t_RFRAC) pari_err(operf,"g",x,y);
2832 p1 = gdeuc(gel(y,2), ggcd(gel(x,2), gel(y,2)));
2833 tetpil = avma;
2834 gel(z,2) = gerepile((pari_sp)z,tetpil,gmul(p1, gel(x,2)));
2835 gel(z,1) = ggcd(gel(x,1), gel(y,1)); return z;
2836 }
2837 pari_err(operf,"g",x,y);
2838 return NULL; /* not reached */
2839 }
2840
2841 /* x a t_VEC,t_COL or t_MAT */
2842 static GEN
vec_lcm(GEN x)2843 vec_lcm(GEN x)
2844 {
2845 if (typ(x) == t_MAT)
2846 {
2847 long i, l = lg(x);
2848 GEN z = cgetg(l, t_VEC);
2849 for (i = 1; i < l; i++) gel(z,i) = glcm0(gel(x,i), NULL);
2850 x = z;
2851 }
2852 return glcm0(x, NULL);
2853 }
2854 static GEN
scal_lcm(GEN x,GEN y)2855 scal_lcm(GEN x, GEN y)
2856 {
2857 pari_sp av = avma;
2858 long tx = typ(x), ty = typ(y);
2859 if (is_matvec_t(tx)) x = vec_lcm(x);
2860 if (is_matvec_t(ty)) y = vec_lcm(y);
2861 return gerepileupto(av, glcm(x, y));
2862 }
2863
2864 static GEN
fix_lcm(GEN x)2865 fix_lcm(GEN x)
2866 {
2867 GEN t;
2868 switch(typ(x))
2869 {
2870 case t_INT: if (signe(x)<0) x = negi(x);
2871 break;
2872 case t_POL:
2873 if (lg(x) <= 2) break;
2874 t = leading_term(x);
2875 if (typ(t) == t_INT && signe(t) < 0) x = gneg(x);
2876 }
2877 return x;
2878 }
2879
2880 GEN
glcm0(GEN x,GEN y)2881 glcm0(GEN x, GEN y) {
2882 if (!y && lg(x) == 2)
2883 {
2884 long tx = typ(x);
2885 if (is_vec_t(tx))
2886 {
2887 x = gel(x,1);
2888 tx = typ(x);
2889 return is_matvec_t(tx)? vec_lcm(x): fix_lcm(x);
2890 }
2891 }
2892 return gassoc_proto(scal_lcm,x,y);
2893 }
2894
2895 GEN
glcm(GEN x,GEN y)2896 glcm(GEN x, GEN y)
2897 {
2898 long tx, ty, i, l;
2899 pari_sp av;
2900 GEN p1,p2,z;
2901
2902 ty = typ(y);
2903 if (is_matvec_t(ty))
2904 {
2905 l=lg(y); z=cgetg(l,ty);
2906 for (i=1; i<l; i++) gel(z,i) = glcm(x,gel(y,i));
2907 return z;
2908 }
2909 tx = typ(x);
2910 if (is_matvec_t(tx))
2911 {
2912 l=lg(x); z=cgetg(l,tx);
2913 for (i=1; i<l; i++) gel(z,i) = glcm(gel(x,i),y);
2914 return z;
2915 }
2916 if (tx==t_INT && ty==t_INT) return lcmii(x,y);
2917 if (gcmp0(x)) return gen_0;
2918
2919 av = avma;
2920 p1 = ggcd(x,y); if (!gcmp1(p1)) y = gdiv(y,p1);
2921 p2 = gmul(x,y);
2922 return gerepileupto(av,fix_lcm(p2));
2923 }
2924
2925 /* x + r ~ x ? Assume x,r are t_POL, deg(r) <= deg(x) */
2926 static int
pol_approx0(GEN r,GEN x,int exact)2927 pol_approx0(GEN r, GEN x, int exact)
2928 {
2929 long i, lx,lr;
2930 if (exact) return gcmp0(r);
2931 lx = lg(x);
2932 lr = lg(r); if (lr < lx) lx = lr;
2933 for (i=2; i<lx; i++)
2934 if (!approx_0(gel(r,i), gel(x,i))) return 0;
2935 return 1;
2936 }
2937
2938 GEN
RgX_gcd_simple(GEN x,GEN y)2939 RgX_gcd_simple(GEN x, GEN y)
2940 {
2941 pari_sp av1, av = avma, lim = stack_lim(av, 1);
2942 GEN r, yorig = y;
2943 int exact = !(isinexactreal(x) || isinexactreal(y));
2944
2945 for(;;)
2946 {
2947 av1 = avma; r = grem(x,y);
2948 if (pol_approx0(r, x, exact))
2949 {
2950 avma = av1;
2951 if (y == yorig) return gcopy(y);
2952 y = normalizepol_approx(y, lg(y));
2953 if (lg(y) == 3) { avma = av; return gen_1; }
2954 return gerepileupto(av,y);
2955 }
2956 x = y; y = r;
2957 if (low_stack(lim,stack_lim(av,1))) {
2958 if(DEBUGMEM>1) pari_warn(warnmem,"RgX_gcd_simple");
2959 gerepileall(av,2, &x,&y);
2960 }
2961 }
2962 }
2963 GEN
RgX_extgcd_simple(GEN a,GEN b,GEN * pu,GEN * pv)2964 RgX_extgcd_simple(GEN a, GEN b, GEN *pu, GEN *pv)
2965 {
2966 pari_sp av = avma;
2967 GEN q, r, d, d1, u, v, v1;
2968 int exact = !(isinexactreal(a) || isinexactreal(b));
2969
2970 d = a; d1 = b; v = gen_0; v1 = gen_1;
2971 for(;;)
2972 {
2973 if (pol_approx0(d1, a, exact)) break;
2974 q = poldivrem(d,d1, &r);
2975 v = gadd(v, gneg_i(gmul(q,v1)));
2976 u=v; v=v1; v1=u;
2977 u=r; d=d1; d1=u;
2978 }
2979 u = gadd(d, gneg_i(gmul(b,v)));
2980 u = RgX_div(u,a);
2981
2982 gerepileall(av, 3, &u,&v,&d);
2983 return d;
2984 }
2985
2986 static int issimplefield(GEN x);
2987 static int
issimplepol(GEN x)2988 issimplepol(GEN x)
2989 {
2990 long i, lx = lg(x);
2991 for (i=2; i<lx; i++)
2992 if (issimplefield(gel(x,i))) return 1;
2993 return 0;
2994 }
2995 /* return 1 if coeff explosion is not possible */
2996 static int
issimplefield(GEN x)2997 issimplefield(GEN x)
2998 {
2999 switch(typ(x))
3000 {
3001 case t_REAL: case t_INTMOD: case t_PADIC: case t_SER:
3002 return 1;
3003 case t_COMPLEX:
3004 return issimplefield(gel(x,1)) || issimplefield(gel(x,2));
3005 case t_POLMOD:
3006 return (typ(x[2]) == t_POL && issimplepol(gel(x,2)))
3007 || issimplefield(gel(x,2)) || issimplepol(gel(x,1));
3008 }
3009 return 0;
3010 }
3011
3012 static int
can_use_modular_gcd(GEN x)3013 can_use_modular_gcd(GEN x)
3014 {
3015 long i;
3016 for (i = lg(x)-1; i > 1; i--)
3017 {
3018 long t = typ(gel(x,i));
3019 if (!is_rational_t(t)) return 0;
3020 }
3021 return 1;
3022 }
3023
3024 static GEN
gcdmonome(GEN x,GEN y)3025 gcdmonome(GEN x, GEN y)
3026 {
3027 long dx=degpol(x), v=varn(x), e=gval(y, v);
3028 pari_sp av = avma;
3029 GEN t = ggcd(leading_term(x), content(y));
3030
3031 if (dx < e) e = dx;
3032 return gerepileupto(av, monomialcopy(t, e, v));
3033 }
3034
3035 /*******************************************************************/
3036 /* */
3037 /* CONTENT / PRIMITIVE PART */
3038 /* */
3039 /*******************************************************************/
3040
3041 GEN
content(GEN x)3042 content(GEN x)
3043 {
3044 long lx, i, t, tx = typ(x);
3045 pari_sp av = avma;
3046 GEN c;
3047
3048 if (is_scalar_t(tx)) return zero_gcd(x, tx);
3049 switch(tx)
3050 {
3051 case t_RFRAC:
3052 {
3053 GEN n = gel(x,1), d = gel(x,2);
3054 /* -- varncmp(vn, vd) < 0 can't happen
3055 * -- if n is POLMOD, its main variable (in the sense of gvar2)
3056 * has lower priority than denominator */
3057 if (typ(n) == t_POLMOD || varncmp(gvar(n), varn(d)) > 0)
3058 n = isinexact(n)? zero_gcd(n, typ(n)): gcopy(n);
3059 else
3060 n = content(n);
3061 return gerepileupto(av, gdiv(n, content(d)));
3062 }
3063
3064 case t_VEC: case t_COL:
3065 lx = lg(x); if (lx==1) return gen_1;
3066 break;
3067
3068 case t_MAT:
3069 {
3070 long hx, j;
3071 lx = lg(x);
3072 if (lx == 1) return gen_1;
3073 hx = lg(x[1]);
3074 if (hx == 1) return gen_1;
3075 if (lx == 2) { x = gel(x,1); lx = lg(x); break; }
3076 if (hx == 2) { x = row_i(x, 1, 1, lx-1); break; }
3077 c = content(gel(x,1));
3078 for (j=2; j<lx; j++)
3079 for (i=1; i<hx; i++) c = ggcd(c,gcoeff(x,i,j));
3080 if (typ(c) == t_INTMOD || isinexact(c)) { avma=av; return gen_1; }
3081 return gerepileupto(av,c);
3082 }
3083
3084 case t_POL: case t_SER:
3085 lx = lg(x); if (lx == 2) return gen_0;
3086 break;
3087 case t_QFR: case t_QFI:
3088 lx = 4; break;
3089
3090 default: pari_err(typeer,"content");
3091 return NULL; /* not reached */
3092 }
3093 for (i=lontyp[tx]; i<lx; i++)
3094 if (typ(x[i]) != t_INT) break;
3095 lx--; c = gel(x,lx);
3096 t = typ(c); if (is_matvec_t(t)) c = content(c);
3097 if (i > lx)
3098 { /* integer coeffs */
3099 while (lx-- > lontyp[tx])
3100 {
3101 c = gcdii(c, gel(x,lx));
3102 if (is_pm1(c)) { avma=av; return gen_1; }
3103 }
3104 }
3105 else
3106 {
3107 if (isinexact(c)) c = zero_gcd(c, typ(c));
3108 while (lx-- > lontyp[tx])
3109 {
3110 GEN d = gel(x,lx);
3111 t = typ(d); if (is_matvec_t(t)) d = content(d);
3112 c = ggcd(c, d);
3113 }
3114 if (typ(c) == t_INTMOD || isinexact(c)) { avma=av; return gen_1; }
3115 }
3116 switch(typ(c))
3117 {
3118 case t_INT:
3119 if (signe(c) < 0) c = negi(c);
3120 break;
3121 case t_VEC: case t_COL: case t_MAT:
3122 pari_err(typeer, "content");
3123 }
3124
3125 return av==avma? gcopy(c): gerepileupto(av,c);
3126 }
3127
3128 GEN
primitive_part(GEN x,GEN * ptc)3129 primitive_part(GEN x, GEN *ptc)
3130 {
3131 pari_sp av = avma;
3132 GEN c = content(x);
3133 if (gcmp1(c)) { avma = av; c = NULL; }
3134 else if (!gcmp0(c)) x = gdiv(x,c);
3135 if (ptc) *ptc = c;
3136 return x;
3137 }
3138
3139 /* NOT MEMORY CLEAN
3140 * As content(), but over Q. Treats polynomial as elts of Q[x1,...xn], instead
3141 * of Q(x2,...,xn)[x1] */
3142 GEN
Q_content(GEN x)3143 Q_content(GEN x)
3144 {
3145 long i, l = lg(x);
3146 GEN d;
3147 pari_sp av;
3148
3149 switch(typ(x))
3150 {
3151 case t_INT: return absi(x);
3152 case t_FRAC: return gabs(x,0);
3153
3154 case t_VEC: case t_COL: case t_MAT:
3155 l = lg(x); if (l == 1) return gen_1;
3156 av = avma; d = Q_content(gel(x,1));
3157 for (i=2; i<l; i++)
3158 {
3159 d = ggcd(d, Q_content(gel(x,i)));
3160 if ((i & 255) == 0) d = gerepileupto(av, d);
3161 }
3162 return gerepileupto(av, d);
3163
3164 case t_POL:
3165 l = lg(x); if (l == 2) return gen_0;
3166 av = avma; d = Q_content(gel(x,2));
3167 for (i=3; i<l; i++)
3168 d = ggcd(d, Q_content(gel(x,i)));
3169 return gerepileupto(av, d);
3170 case t_POLMOD: return Q_content(gel(x,2));
3171 case t_COMPLEX: return ggcd(Q_content(gel(x,1)), Q_content(gel(x,2)));
3172 }
3173 pari_err(typeer,"Q_content");
3174 return NULL; /* not reached */
3175 }
3176
3177 GEN
Q_primitive_part(GEN x,GEN * ptc)3178 Q_primitive_part(GEN x, GEN *ptc)
3179 {
3180 pari_sp av = avma;
3181 GEN c = Q_content(x);
3182 if (gcmp1(c)) { avma = av; c = NULL; }
3183 else if (!gcmp0(c)) x = Q_div_to_int(x,c);
3184 if (ptc) *ptc = c;
3185 return x;
3186 }
3187
3188 GEN
primpart(GEN x)3189 primpart(GEN x) { return primitive_part(x, NULL); }
3190
3191 GEN
Q_primpart(GEN x)3192 Q_primpart(GEN x) { return Q_primitive_part(x, NULL); }
3193
3194 /* NOT MEMORY CLEAN (because of t_FRAC).
3195 * As denom(), but over Q. Treats polynomial as elts of Q[x1,...xn], instead
3196 * of Q(x2,...,xn)[x1] */
3197 GEN
Q_denom(GEN x)3198 Q_denom(GEN x)
3199 {
3200 long i, l = lg(x);
3201 GEN d, D;
3202 pari_sp av;
3203
3204 switch(typ(x))
3205 {
3206 case t_INT: return gen_1;
3207 case t_FRAC: return gel(x,2);
3208
3209 case t_VEC: case t_COL: case t_MAT:
3210 l = lg(x); if (l == 1) return gen_1;
3211 av = avma; d = Q_denom(gel(x,1));
3212 for (i=2; i<l; i++)
3213 {
3214 D = Q_denom(gel(x,i));
3215 if (D != gen_1) d = lcmii(d, D);
3216 if ((i & 255) == 0) d = gerepileuptoint(av, d);
3217 }
3218 return gerepileuptoint(av, d);
3219
3220 case t_POL:
3221 l = lg(x); if (l == 2) return gen_1;
3222 av = avma; d = Q_denom(gel(x,2));
3223 for (i=3; i<l; i++)
3224 {
3225 D = Q_denom(gel(x,i));
3226 if (D != gen_1) d = lcmii(d, D);
3227 }
3228 return gerepileuptoint(av, d);
3229 }
3230 pari_err(typeer,"Q_denom");
3231 return NULL; /* not reached */
3232 }
3233
3234 GEN
Q_remove_denom(GEN x,GEN * ptd)3235 Q_remove_denom(GEN x, GEN *ptd)
3236 {
3237 GEN d = Q_denom(x);
3238 if (gcmp1(d)) d = NULL; else x = Q_muli_to_int(x,d);
3239 if (ptd) *ptd = d;
3240 return x;
3241 }
3242
3243 /* return y = x * d, assuming x rational, and d,y integral */
3244 GEN
Q_muli_to_int(GEN x,GEN d)3245 Q_muli_to_int(GEN x, GEN d)
3246 {
3247 long i, l, t = typ(x);
3248 GEN y, xn, xd;
3249 pari_sp av;
3250
3251 if (typ(d) != t_INT) pari_err(typeer,"Q_muli_to_int");
3252 switch (t)
3253 {
3254 case t_INT:
3255 return mulii(x,d);
3256
3257 case t_FRAC:
3258 xn = gel(x,1);
3259 xd = gel(x,2); av = avma;
3260 y = mulii(xn, diviiexact(d, xd));
3261 return gerepileuptoint(av, y);
3262
3263 case t_VEC: case t_COL: case t_MAT:
3264 l = lg(x); y = cgetg(l, t);
3265 for (i=1; i<l; i++) gel(y,i) = Q_muli_to_int(gel(x,i), d);
3266 return y;
3267
3268 case t_POL:
3269 l = lg(x); y = cgetg(l, t_POL); y[1] = x[1];
3270 for (i=2; i<l; i++) gel(y,i) = Q_muli_to_int(gel(x,i), d);
3271 return y;
3272
3273 case t_POLMOD:
3274 y = cgetg(3, t_POLMOD);
3275 gel(y,1) = gcopy(gel(x,1));
3276 gel(y,2) = Q_muli_to_int(gel(x,2), d);
3277 return y;
3278 }
3279 pari_err(typeer,"Q_muli_to_int");
3280 return NULL; /* not reached */
3281 }
3282
3283 /* return x * n/d. x: rational; d,n,result: integral. n = NULL represents 1 */
3284 GEN
Q_divmuli_to_int(GEN x,GEN d,GEN n)3285 Q_divmuli_to_int(GEN x, GEN d, GEN n)
3286 {
3287 long i, l, t = typ(x);
3288 GEN y, xn, xd;
3289 pari_sp av;
3290
3291 switch(t)
3292 {
3293 case t_INT:
3294 av = avma; y = diviiexact(x,d);
3295 if (n) y = gerepileuptoint(av, mulii(y,n));
3296 return y;
3297
3298 case t_FRAC:
3299 xn = gel(x,1);
3300 xd = gel(x,2); av = avma;
3301 y = mulii(diviiexact(xn, d), diviiexact(n, xd));
3302 return gerepileuptoint(av, y);
3303
3304 case t_VEC: case t_COL: case t_MAT:
3305 l = lg(x); y = cgetg(l, t);
3306 for (i=1; i<l; i++) gel(y,i) = Q_divmuli_to_int(gel(x,i), d,n);
3307 return y;
3308
3309 case t_POL:
3310 l = lg(x); y = cgetg(l, t_POL); y[1] = x[1];
3311 for (i=2; i<l; i++) gel(y,i) = Q_divmuli_to_int(gel(x,i), d,n);
3312 return y;
3313
3314 case t_POLMOD:
3315 y = cgetg(3, t_POLMOD);
3316 gel(y,1) = gcopy(gel(x,1));
3317 gel(y,2) = Q_divmuli_to_int(gel(x,2), d,n);
3318 return y;
3319 }
3320 pari_err(typeer,"Q_divmuli_to_int");
3321 return NULL; /* not reached */
3322 }
3323
3324 /* return y = x / c, assuming x,c rational, and y integral */
3325 GEN
Q_div_to_int(GEN x,GEN c)3326 Q_div_to_int(GEN x, GEN c)
3327 {
3328 GEN d, n;
3329 switch(typ(c))
3330 {
3331 case t_INT:
3332 return Q_divmuli_to_int(x, c, NULL);
3333 case t_FRAC:
3334 n = gel(c,1);
3335 d = gel(c,2); if (gcmp1(n)) return Q_muli_to_int(x,d);
3336 return Q_divmuli_to_int(x, n,d);
3337 }
3338 pari_err(typeer,"Q_div_to_int");
3339 return NULL; /* not reached */
3340 }
3341
3342 /*******************************************************************/
3343 /* */
3344 /* SUBRESULTANT */
3345 /* */
3346 /*******************************************************************/
3347 /* for internal use */
3348 GEN
gdivexact(GEN x,GEN y)3349 gdivexact(GEN x, GEN y)
3350 {
3351 long i,lx;
3352 GEN z;
3353 if (gcmp1(y)) return x;
3354 switch(typ(x))
3355 {
3356 case t_INT:
3357 if (typ(y)==t_INT) return diviiexact(x,y);
3358 if (!signe(x)) return gen_0;
3359 break;
3360 case t_INTMOD:
3361 case t_POLMOD: return gmul(x,ginv(y));
3362 case t_POL:
3363 switch(typ(y))
3364 {
3365 case t_INTMOD:
3366 case t_POLMOD: return gmul(x,ginv(y));
3367 case t_POL:
3368 if (varn(x)==varn(y)) return poldivrem(x,y, NULL);
3369 }
3370 lx = lg(x); z = new_chunk(lx);
3371 for (i=2; i<lx; i++) gel(z,i) = gdivexact(gel(x,i),y);
3372 z[1] = x[1];
3373 z[0] = x[0]; return z;
3374 case t_VEC: case t_COL: case t_MAT:
3375 lx = lg(x); z = new_chunk(lx);
3376 for (i=1; i<lx; i++) gel(z,i) = gdivexact(gel(x,i),y);
3377 z[0] = x[0]; return z;
3378 }
3379 if (DEBUGLEVEL) pari_warn(warner,"missing case in gdivexact");
3380 return gdiv(x,y);
3381 }
3382
3383 static GEN
init_resultant(GEN x,GEN y)3384 init_resultant(GEN x, GEN y)
3385 {
3386 long tx,ty;
3387 if (gcmp0(x) || gcmp0(y)) return gen_0;
3388 tx = typ(x); ty = typ(y);
3389 if (is_scalar_t(tx) || is_scalar_t(ty))
3390 {
3391 if (tx==t_POL) return gpowgs(y,degpol(x));
3392 if (ty==t_POL) return gpowgs(x,degpol(y));
3393 return gen_1;
3394 }
3395 if (tx!=t_POL || ty!=t_POL) pari_err(typeer,"subresall");
3396 if (varn(x)==varn(y)) return NULL;
3397 return (varn(x)<varn(y))? gpowgs(y,degpol(x)): gpowgs(x,degpol(y));
3398 }
3399
3400 /* return coefficients s.t x = x_0 X^n + ... + x_n */
3401 GEN
revpol(GEN x)3402 revpol(GEN x)
3403 {
3404 long i,n = degpol(x);
3405 GEN y = cgetg(n+3,t_POL);
3406 y[1] = x[1]; x += 2; y += 2;
3407 for (i=0; i<=n; i++) y[i] = x[n-i];
3408 return y;
3409 }
3410
3411 /* assume dx >= dy, y non constant, mod either NULL or a t_POL. */
3412 static GEN
pseudorem_i(GEN x,GEN y,GEN mod)3413 pseudorem_i(GEN x, GEN y, GEN mod)
3414 {
3415 long vx = varn(x), dx, dy, dz, i, lx, p;
3416 pari_sp av = avma, av2, lim;
3417
3418 if (!signe(y)) pari_err(gdiver);
3419 (void)new_chunk(2);
3420 dx=degpol(x); x = revpol(x);
3421 dy=degpol(y); y = revpol(y); dz=dx-dy; p = dz+1;
3422 av2 = avma; lim = stack_lim(av2,1);
3423 for (;;)
3424 {
3425 gel(x,0) = gneg(gel(x,0)); p--;
3426 for (i=1; i<=dy; i++)
3427 {
3428 gel(x,i) = gadd(gmul(gel(y,0), gel(x,i)), gmul(gel(x,0),gel(y,i)));
3429 if (mod) gel(x,i) = RgX_rem(gel(x,i), mod);
3430 }
3431 for ( ; i<=dx; i++)
3432 {
3433 gel(x,i) = gmul(gel(y,0), gel(x,i));
3434 if (mod) gel(x,i) = RgX_rem(gel(x,i), mod);
3435 }
3436 do { x++; dx--; } while (dx >= 0 && gcmp0(gel(x,0)));
3437 if (dx < dy) break;
3438 if (low_stack(lim,stack_lim(av2,1)))
3439 {
3440 if(DEBUGMEM>1) pari_warn(warnmem,"pseudorem dx = %ld >= %ld",dx,dy);
3441 gerepilecoeffs(av2,x,dx+1);
3442 }
3443 }
3444 if (dx < 0) return zeropol(vx);
3445 lx = dx+3; x -= 2;
3446 x[0] = evaltyp(t_POL) | evallg(lx);
3447 x[1] = evalsigne(1) | evalvarn(vx);
3448 x = revpol(x) - 2;
3449 if (p)
3450 { /* multiply by y[0]^p [beware dummy vars from FpY_FpXY_resultant] */
3451 GEN t = gel(y,0);
3452 if (mod)
3453 { /* assume p fairly small */
3454 for (i=1; i<p; i++)
3455 t = RgX_rem(gmul(t, gel(y,0)), mod);
3456 }
3457 else
3458 t = gpowgs(t, p);
3459 for (i=2; i<lx; i++)
3460 {
3461 gel(x,i) = gmul(gel(x,i), t);
3462 if (mod) gel(x,i) = RgX_rem(gel(x,i), mod);
3463 }
3464 if (!mod) return gerepileupto(av, x);
3465 }
3466 return gerepilecopy(av, x);
3467 }
3468
3469 GEN
pseudorem(GEN x,GEN y)3470 pseudorem(GEN x, GEN y) { return pseudorem_i(x,y, NULL); }
3471
3472 /* assume dx >= dy, y non constant
3473 * Compute z,r s.t lc(y)^(dx-dy+1) x = z y + r */
3474 GEN
pseudodiv(GEN x,GEN y,GEN * ptr)3475 pseudodiv(GEN x, GEN y, GEN *ptr)
3476 {
3477 long vx = varn(x), dx, dy, dz, i, iz, lx, lz, p;
3478 pari_sp av = avma, av2, lim;
3479 GEN z, r, ypow;
3480
3481 if (!signe(y)) pari_err(gdiver);
3482 (void)new_chunk(2);
3483 dx=degpol(x); x = revpol(x);
3484 dy=degpol(y); y = revpol(y); dz=dx-dy; p = dz+1;
3485 lz = dz+3; z = cgetg(lz, t_POL) + 2;
3486 ypow = new_chunk(dz+1);
3487 gel(ypow,0) = gen_1;
3488 for (i=1; i<=dz; i++) gel(ypow,i) = gmul(gel(ypow,i-1), gel(y,0));
3489 av2 = avma; lim = stack_lim(av2,1);
3490 for (iz=0;;)
3491 {
3492 p--;
3493 gel(z,iz++) = gmul(gel(x,0), gel(ypow,p));
3494 gel(x,0) = gneg(gel(x,0));
3495 for (i=1; i<=dy; i++)
3496 gel(x,i) = gadd(gmul(gel(y,0), gel(x,i)), gmul(gel(x,0),gel(y,i)));
3497 for ( ; i<=dx; i++)
3498 gel(x,i) = gmul(gel(y,0), gel(x,i));
3499 x++; dx--;
3500 while (dx >= dy && gcmp0(gel(x,0))) { x++; dx--; gel(z,iz++) = gen_0; }
3501 if (dx < dy) break;
3502 if (low_stack(lim,stack_lim(av2,1)))
3503 {
3504 if(DEBUGMEM>1) pari_warn(warnmem,"pseudodiv dx = %ld >= %ld",dx,dy);
3505 gerepilecoeffs2(av2,x,dx+1, z,iz);
3506 }
3507 }
3508 while (dx >= 0 && gcmp0(gel(x,0))) { x++; dx--; }
3509 if (dx < 0)
3510 x = zeropol(vx);
3511 else
3512 {
3513 lx = dx+3; x -= 2;
3514 x[0] = evaltyp(t_POL) | evallg(lx);
3515 x[1] = evalsigne(1) | evalvarn(vx);
3516 x = revpol(x) - 2;
3517 }
3518
3519 z -= 2;
3520 z[0] = evaltyp(t_POL) | evallg(lz);
3521 z[1] = evalsigne(1) | evalvarn(vx);
3522 z = revpol(z) - 2;
3523 r = gmul(x, gel(ypow,p));
3524 gerepileall(av, 2, &z, &r);
3525 *ptr = r; return z;
3526 }
3527
3528 /* Return resultant(u,v). If sol != NULL: set *sol to the last non-zero
3529 * polynomial in the prs IF the sequence was computed, and gen_0 otherwise */
3530 GEN
subresall(GEN u,GEN v,GEN * sol)3531 subresall(GEN u, GEN v, GEN *sol)
3532 {
3533 pari_sp av, av2, lim;
3534 long degq,dx,dy,du,dv,dr,signh;
3535 GEN z,g,h,r,p1,p2,cu,cv;
3536
3537 if (sol) *sol=gen_0;
3538 if ((r = init_resultant(u,v))) return r;
3539
3540 if (isinexact(u) || isinexact(v)) return resultant2(u,v);
3541 dx=degpol(u); dy=degpol(v); signh=1;
3542 if (dx < dy)
3543 {
3544 swap(u,v); lswap(dx,dy);
3545 if (both_odd(dx, dy)) signh = -signh;
3546 }
3547 if (dy==0) return gpowgs(gel(v,2),dx);
3548 av = avma;
3549 u = primitive_part(u, &cu);
3550 v = primitive_part(v, &cv);
3551 g = h = gen_1; av2 = avma; lim = stack_lim(av2,1);
3552 for(;;)
3553 {
3554 r = pseudorem(u,v); dr = lg(r);
3555 if (dr == 2)
3556 {
3557 if (sol) { avma = (pari_sp)(r+2); *sol=gerepileupto(av,v); } else avma = av;
3558 return gen_0;
3559 }
3560 du = degpol(u); dv = degpol(v); degq = du-dv;
3561 u = v; p1 = g; g = leading_term(u);
3562 switch(degq)
3563 {
3564 case 0: break;
3565 case 1:
3566 p1 = gmul(h,p1); h = g; break;
3567 default:
3568 p1 = gmul(gpowgs(h,degq),p1);
3569 h = gdivexact(gpowgs(g,degq), gpowgs(h,degq-1));
3570 }
3571 if (both_odd(du,dv)) signh = -signh;
3572 v = gdivexact(r,p1);
3573 if (dr==3) break;
3574 if (low_stack(lim,stack_lim(av2,1)))
3575 {
3576 if(DEBUGMEM>1) pari_warn(warnmem,"subresall, dr = %ld",dr);
3577 gerepileall(av2,4, &u, &v, &g, &h);
3578 }
3579 }
3580 z = gel(v,2);
3581 if (dv > 1) z = gdivexact(gpowgs(z,dv), gpowgs(h,dv-1));
3582 if (signh < 0) z = gneg(z); /* z = resultant(ppart(x), ppart(y)) */
3583 p2 = gen_1;
3584 if (cu) p2 = gmul(p2, gpowgs(cu,dy));
3585 if (cv) p2 = gmul(p2, gpowgs(cv,dx));
3586 z = gmul(z, p2);
3587
3588 if (sol) u = gclone(u);
3589 z = gerepileupto(av, z);
3590 if (sol) { *sol = gcopy(u); gunclone(u); }
3591 return z;
3592 }
3593
3594 static GEN
scalar_res(GEN x,GEN y,GEN * U,GEN * V)3595 scalar_res(GEN x, GEN y, GEN *U, GEN *V)
3596 {
3597 *V=gpowgs(y,degpol(x)-1); *U=gen_0; return gmul(y,*V);
3598 }
3599
3600 /* compute U, V s.t Ux + Vy = resultant(x,y) */
3601 GEN
subresext(GEN x,GEN y,GEN * U,GEN * V)3602 subresext(GEN x, GEN y, GEN *U, GEN *V)
3603 {
3604 pari_sp av, av2, tetpil, lim;
3605 long dx, dy, signh, tx = typ(x), ty = typ(y);
3606 GEN z, g, h, p1, cu, cv, u, v, um1, uze, vze;
3607 GEN *gptr[3];
3608
3609 if (!is_extscalar_t(tx) || !is_extscalar_t(ty)) pari_err(typeer,"subresext");
3610 if (gcmp0(x) || gcmp0(y)) { *U = *V = gen_0; return gen_0; }
3611 if (tx != t_POL) {
3612 if (ty != t_POL) { *U = ginv(x); *V = gen_0; return gen_1; }
3613 return scalar_res(y,x,V,U);
3614 }
3615 if (ty != t_POL) return scalar_res(x,y,U,V);
3616 if (varn(x) != varn(y))
3617 return varncmp(varn(x), varn(y)) < 0? scalar_res(x,y,U,V)
3618 : scalar_res(y,x,V,U);
3619 dx = degpol(x); dy = degpol(y); signh = 1;
3620 if (dx < dy)
3621 {
3622 pswap(U,V); lswap(dx,dy); swap(x,y);
3623 if (both_odd(dx, dy)) signh = -signh;
3624 }
3625 if (dy == 0)
3626 {
3627 *V = gpowgs(gel(y,2),dx-1);
3628 *U = gen_0; return gmul(*V,gel(y,2));
3629 }
3630 av = avma;
3631 u = x = primitive_part(x, &cu);
3632 v = y = primitive_part(y, &cv);
3633 g = h = gen_1; av2 = avma; lim = stack_lim(av2,1);
3634 um1 = gen_1; uze = gen_0;
3635 for(;;)
3636 {
3637 GEN r, q = pseudodiv(u,v, &r);
3638 long du, dv, degq, dr = lg(r);
3639 if (dr == 2) { *U = *V = gen_0; avma = av; return gen_0; }
3640
3641 du = degpol(u); dv = degpol(v); degq = du-dv;
3642 /* lead(v)^(degq + 1) * um1 - q * uze */
3643 p1 = gsub(gmul(gpowgs(gel(v,dv+2),degq+1),um1), gmul(q,uze));
3644 um1 = uze; uze = p1;
3645 u = v; p1 = g; g = leading_term(u);
3646 switch(degq)
3647 {
3648 case 0: break;
3649 case 1: p1 = gmul(h,p1); h = g; break;
3650 default:
3651 p1 = gmul(gpowgs(h,degq),p1);
3652 h = gdivexact(gpowgs(g,degq), gpowgs(h,degq-1));
3653 }
3654 if (both_odd(du, dv)) signh = -signh;
3655 v = gdivexact(r,p1);
3656 uze= gdivexact(uze,p1);
3657 if (dr == 3) {
3658 z = gel(v,2);
3659 if (dv > 1)
3660 { /* z = gdivexact(gpowgs(z,dv), gpowgs(h,dv-1)); */
3661 p1 = gpowgs(gdiv(z,h),dv-1);
3662 z = gmul(z,p1); uze = gmul(uze,p1);
3663 }
3664 break;
3665 }
3666 if (low_stack(lim,stack_lim(av2,1)))
3667 {
3668 if(DEBUGMEM>1) pari_warn(warnmem,"subresext, dr = %ld",dr);
3669 gerepileall(av2,6, &u,&v,&g,&h,&uze,&um1);
3670 }
3671 }
3672 if (signh < 0) { z = gneg_i(z); uze = gneg_i(uze); }
3673 p1 = gadd(z, gneg(gmul(uze,x)));
3674 vze = RgX_divrem(p1, y, &p1);
3675 if (!gcmp0(p1)) pari_warn(warner,"inexact computation in subresext");
3676 /* uze ppart(x) + vze ppart(y) = z = resultant(ppart(x), ppart(y)), */
3677 p1 = gen_1;
3678 if (cu) p1 = gmul(p1, gpowgs(cu,dy));
3679 if (cv) p1 = gmul(p1, gpowgs(cv,dx));
3680 cu = cu? gdiv(p1,cu): p1;
3681 cv = cv? gdiv(p1,cv): p1;
3682
3683 tetpil = avma;
3684 z = gmul(z,p1);
3685 *U = gmul(uze,cu);
3686 *V = gmul(vze,cv);
3687 gptr[0]=&z; gptr[1]=U; gptr[2]=V;
3688 gerepilemanysp(av,tetpil,gptr,3);
3689 return z;
3690 }
3691
3692 static GEN
scalar_bezout(GEN x,GEN y,GEN * U,GEN * V)3693 scalar_bezout(GEN x, GEN y, GEN *U, GEN *V)
3694 {
3695 *U=gen_0; *V=ginv(y); return pol_1[varn(x)];
3696 }
3697
3698 static GEN
zero_bezout(GEN y,GEN * U,GEN * V)3699 zero_bezout(GEN y, GEN *U, GEN *V)
3700 {
3701 GEN x=content(y);
3702 *U=gen_0; *V = ginv(x); return gmul(y,*V);
3703 }
3704
3705 /* compute U, V s.t Ux + Vy = GCD(x,y) using subresultant */
3706 GEN
RgX_extgcd(GEN x,GEN y,GEN * U,GEN * V)3707 RgX_extgcd(GEN x, GEN y, GEN *U, GEN *V)
3708 {
3709 pari_sp av, av2, tetpil, lim;
3710 long dx, dy, tx = typ(x), ty = typ(y);
3711 GEN z, g, h, p1, cu, cv, u, v, um1, uze, vze, *gptr[3];
3712
3713 if (!is_extscalar_t(tx) || !is_extscalar_t(ty)) pari_err(typeer,"subresext");
3714 if (gcmp0(x)) {
3715 if (gcmp0(y)) { *U = *V = gen_0; return gen_0; }
3716 return zero_bezout(y,U,V);
3717 }
3718 if (gcmp0(y)) return zero_bezout(x,V,U);
3719 if (tx != t_POL) {
3720 if (ty != t_POL) { *U = ginv(x); *V = gen_0; return pol_1[0]; }
3721 return scalar_bezout(y,x,V,U);
3722 }
3723 if (ty != t_POL) return scalar_bezout(x,y,U,V);
3724 if (varn(x) != varn(y))
3725 return varncmp(varn(x), varn(y)) < 0? scalar_bezout(x,y,U,V)
3726 : scalar_bezout(y,x,V,U);
3727 dx = degpol(x); dy = degpol(y);
3728 if (dx < dy)
3729 {
3730 pswap(U,V); lswap(dx,dy); swap(x,y);
3731 }
3732 if (dy==0) return scalar_bezout(x,y,U,V);
3733
3734 av = avma;
3735 u = x = primitive_part(x, &cu);
3736 v = y = primitive_part(y, &cv);
3737 g = h = gen_1; av2 = avma; lim = stack_lim(av2,1);
3738 um1 = gen_1; uze = gen_0;
3739 for(;;)
3740 {
3741 GEN r, q = pseudodiv(u,v, &r);
3742 long du, dv, degq, dr = lg(r);
3743 if (dr == 2) break;
3744
3745 du = degpol(u); dv = degpol(v); degq = du-dv;
3746 p1 = gsub(gmul(gpowgs(gel(v,dv+2),degq+1),um1), gmul(q,uze));
3747 um1 = uze; uze = p1;
3748 u = v; p1 = g; g = leading_term(u);
3749 switch(degq)
3750 {
3751 case 0: break;
3752 case 1:
3753 p1 = gmul(h,p1); h = g; break;
3754 default:
3755 p1 = gmul(gpowgs(h,degq), p1);
3756 h = gdiv(gpowgs(g,degq), gpowgs(h,degq-1));
3757 }
3758 v = gdivexact(r,p1);
3759 uze= gdivexact(uze,p1);
3760 if (dr==3) break;
3761 if (low_stack(lim,stack_lim(av2,1)))
3762 {
3763 if(DEBUGMEM>1) pari_warn(warnmem,"RgX_extgcd, dr = %ld",dr);
3764 gerepileall(av2,6,&u,&v,&g,&h,&uze,&um1);
3765 }
3766 }
3767 p1 = gadd(v, gneg(gmul(uze,x)));
3768 vze = RgX_divrem(p1, y, &p1);
3769 if (!gcmp0(p1)) pari_warn(warner,"inexact computation in RgX_extgcd");
3770 if (cu) uze = gdiv(uze,cu);
3771 if (cv) vze = gdiv(vze,cv);
3772 p1 = ginv(content(v));
3773
3774 tetpil = avma;
3775 *U = gmul(uze,p1);
3776 *V = gmul(vze,p1);
3777 z = gmul(v,p1);
3778 gptr[0]=U; gptr[1]=V; gptr[2]=&z;
3779 gerepilemanysp(av,tetpil,gptr,3); return z;
3780 }
3781
3782 /*******************************************************************/
3783 /* */
3784 /* RESULTANT USING DUCOS VARIANT */
3785 /* */
3786 /*******************************************************************/
3787
3788 static GEN
reductum(GEN P)3789 reductum(GEN P)
3790 {
3791 if (signe(P)==0) return P;
3792 return normalizepol_i(shallowcopy(P),lg(P)-1);
3793 }
3794
3795 static GEN
Lazard(GEN x,GEN y,long n)3796 Lazard(GEN x, GEN y, long n)
3797 {
3798 long a, b;
3799 GEN c;
3800
3801 if (n<=1) return x;
3802 a=1; while (n >= (b=a+a)) a=b;
3803 c=x; n-=a;
3804 while (a>1)
3805 {
3806 a>>=1; c=gdivexact(gsqr(c),y);
3807 if (n>=a) { c=gdivexact(gmul(c,x),y); n -= a; }
3808 }
3809 return c;
3810 }
3811
3812 static GEN
Lazard2(GEN F,GEN x,GEN y,long n)3813 Lazard2(GEN F, GEN x, GEN y, long n)
3814 {
3815 if (n<=1) return F;
3816 return gdivexact(gmul(Lazard(x,y,n-1), F), y);
3817 }
3818
3819 /* deg(P) > deg(Q) */
3820 static GEN
nextSousResultant(GEN P,GEN Q,GEN Z,GEN s)3821 nextSousResultant(GEN P, GEN Q, GEN Z, GEN s)
3822 {
3823 GEN p0, q0, z0, H, A;
3824 long pr, p, q, j, v = varn(P);
3825 pari_sp av, lim;
3826
3827 z0 = leading_term(Z);
3828 p = degpol(P); p0 = leading_term(P); P = reductum(P);
3829 q = degpol(Q); q0 = leading_term(Q); Q = reductum(Q);
3830
3831 av = avma; lim = stack_lim(av,1);
3832 H = gneg(reductum(Z));
3833 pr = degpol(P);
3834 A = (q <= pr)? gmul(gel(P,q+2),H): gen_0;
3835 for (j = q+1; j < p; j++)
3836 {
3837 H = (degpol(H) == q-1) ?
3838 addshift(reductum(H), gdivexact(gmul(gneg(gel(H,q+1)),Q), q0)) :
3839 addshift(H, zeropol(v));
3840 if (j <= pr) A = gadd(A,gmul(gel(P,j+2),H));
3841 if (low_stack(lim,stack_lim(av,1)))
3842 {
3843 if(DEBUGMEM>1) pari_warn(warnmem,"nextSousResultant j = %ld/%ld",j,p);
3844 gerepileall(av,2,&A,&H);
3845 }
3846 }
3847 P = normalizepol_i(P, min(pr+3,q+2));
3848 A = gdivexact(gadd(A,gmul(z0,P)), p0);
3849 A = (degpol(H) == q-1) ?
3850 gadd(gmul(q0,addshift(reductum(H),A)), gmul(gneg(gel(H,q+1)),Q)) :
3851 gmul(q0, addshift(H,A));
3852 return gdivexact(A, ((p-q)&1)? s: gneg(s));
3853 }
3854
3855 GEN
resultantducos(GEN P,GEN Q)3856 resultantducos(GEN P, GEN Q)
3857 {
3858 pari_sp av = avma, av2, lim;
3859 long dP,dQ,delta;
3860 GEN cP,cQ,Z,s;
3861
3862 if ((Z = init_resultant(P,Q))) return Z;
3863 dP = degpol(P);
3864 dQ = degpol(Q);
3865 P = primitive_part(P, &cP);
3866 Q = primitive_part(Q, &cQ);
3867 delta = dP - dQ;
3868 if (delta < 0)
3869 {
3870 Z = (dP & dQ & 1)? gneg(Q): Q;
3871 Q = P; P = Z; delta = -delta;
3872 }
3873 s = gen_1;
3874 if (degpol(Q) > 0)
3875 {
3876 av2 = avma; lim = stack_lim(av2,1);
3877 s = gpowgs(leading_term(Q),delta);
3878 Z = Q;
3879 Q = pseudorem(P, gneg(Q));
3880 P = Z;
3881 while(degpol(Q) > 0)
3882 {
3883 if (low_stack(lim,stack_lim(av,1)))
3884 {
3885 if(DEBUGMEM>1) pari_warn(warnmem,"resultantducos, degpol Q = %ld",degpol(Q));
3886 gerepileall(av2,2,&P,&Q); s = leading_term(P);
3887 }
3888 delta = degpol(P) - degpol(Q);
3889 Z = Lazard2(Q, leading_term(Q), s, delta);
3890 Q = nextSousResultant(P, Q, Z, s);
3891 P = Z;
3892 s = leading_term(P);
3893 }
3894 }
3895 if (!signe(Q)) { avma = av; return gen_0; }
3896 if (!degpol(P)){ avma = av; return gen_1; }
3897 s = Lazard(leading_term(Q), s, degpol(P));
3898 if (cP) s = gmul(s, gpowgs(cP,dQ));
3899 if (cQ) s = gmul(s, gpowgs(cQ,dP)); else if (!cP) s = gcopy(s);
3900 return gerepileupto(av, s);
3901 }
3902
3903 /*******************************************************************/
3904 /* */
3905 /* RESULTANT USING SYLVESTER MATRIX */
3906 /* */
3907 /*******************************************************************/
3908 static GEN
_zeropol(void)3909 _zeropol(void)
3910 {
3911 GEN x = cgetg(3,t_POL);
3912 x[1] = 0;
3913 gel(x,2) = gen_0; return x;
3914 }
3915
3916 static GEN
sylvester_col(GEN x,long j,long d,long D)3917 sylvester_col(GEN x, long j, long d, long D)
3918 {
3919 GEN c = cgetg(d+1,t_COL);
3920 long i;
3921 for (i=1; i< j; i++) gel(c,i) = gen_0;
3922 for ( ; i<=D; i++) c[i]=x[D-i+2];
3923 for ( ; i<=d; i++) gel(c,i) = gen_0;
3924 return c;
3925 }
3926
3927 GEN
sylvestermatrix_i(GEN x,GEN y)3928 sylvestermatrix_i(GEN x, GEN y)
3929 {
3930 long j,d,dx,dy;
3931 GEN M;
3932
3933 dx = degpol(x); if (dx < 0) { dx = 0; x = _zeropol(); }
3934 dy = degpol(y); if (dy < 0) { dy = 0; y = _zeropol(); }
3935 d = dx+dy; M = cgetg(d+1,t_MAT);
3936 for (j=1; j<=dy; j++) gel(M,j) = sylvester_col(x,j,d,j+dx);
3937 for (j=1; j<=dx; j++) gel(M,j+dy) = sylvester_col(y,j,d,j+dy);
3938 return M;
3939 }
3940
3941 GEN
sylvestermatrix(GEN x,GEN y)3942 sylvestermatrix(GEN x, GEN y)
3943 {
3944 long i,j,lx;
3945 if (typ(x)!=t_POL || typ(y)!=t_POL) pari_err(typeer,"sylvestermatrix");
3946 if (varn(x) != varn(y))
3947 pari_err(talker,"not the same variables in sylvestermatrix");
3948 x = sylvestermatrix_i(x,y); lx = lg(x);
3949 for (i=1; i<lx; i++)
3950 for (j=1; j<lx; j++) gcoeff(x,i,j) = gcopy(gcoeff(x,i,j));
3951 return x;
3952 }
3953
3954 GEN
resultant2(GEN x,GEN y)3955 resultant2(GEN x, GEN y)
3956 {
3957 pari_sp av;
3958 GEN r;
3959 if ((r = init_resultant(x,y))) return r;
3960 av = avma; return gerepileupto(av,det(sylvestermatrix_i(x,y)));
3961 }
3962
3963 static GEN
fix_pol(GEN x,long v,long * mx)3964 fix_pol(GEN x, long v, long *mx)
3965 {
3966 long vx;
3967 GEN p1;
3968
3969 if (typ(x) == t_POL)
3970 {
3971 vx = varn(x);
3972 if (vx)
3973 {
3974 if (v>=vx) return gsubst(x,v,pol_x[0]);
3975 p1 = cgetg(3,t_POL);
3976 p1[1] = evalvarn(0)|evalsigne(signe(x));
3977 gel(p1,2) = x; return p1;
3978 }
3979 if (v)
3980 {
3981 *mx = 1;
3982 return gsubst(gsubst(x,0,pol_x[MAXVARN]),v,pol_x[0]);
3983 }
3984 }
3985 return x;
3986 }
3987
3988 /* resultant of x and y with respect to variable v, or with respect to their
3989 * main variable if v < 0. */
3990 GEN
polresultant0(GEN x,GEN y,long v,long flag)3991 polresultant0(GEN x, GEN y, long v, long flag)
3992 {
3993 long m = 0;
3994 pari_sp av = avma;
3995
3996 if (v >= 0)
3997 {
3998 x = fix_pol(x,v, &m);
3999 y = fix_pol(y,v, &m);
4000 }
4001 switch(flag)
4002 {
4003 case 0: x=subresall(x,y,NULL); break;
4004 case 1: x=resultant2(x,y); break;
4005 case 2: x=resultantducos(x,y); break;
4006 default: pari_err(flagerr,"polresultant");
4007 }
4008 if (m) x = gsubst(x,MAXVARN,pol_x[0]);
4009 return gerepileupto(av,x);
4010 }
4011
4012 /*******************************************************************/
4013 /* */
4014 /* GCD USING SUBRESULTANT */
4015 /* */
4016 /*******************************************************************/
4017 GEN
srgcd(GEN x,GEN y)4018 srgcd(GEN x, GEN y)
4019 {
4020 long tx = typ(x), ty = typ(y), dx, dy, vx;
4021 pari_sp av, av1, tetpil, lim;
4022 GEN d, g, h, p1, p2, u, v;
4023
4024 if (!signe(y)) return gcopy(x);
4025 if (!signe(x)) return gcopy(y);
4026 if (is_scalar_t(tx) || is_scalar_t(ty)) return gen_1;
4027 if (tx!=t_POL || ty!=t_POL) pari_err(typeer,"srgcd");
4028 vx=varn(x); if (vx!=varn(y)) return gen_1;
4029 if (ismonome(x)) return gcdmonome(x,y);
4030 if (ismonome(y)) return gcdmonome(y,x);
4031 av = avma;
4032 if (can_use_modular_gcd(x) &&
4033 can_use_modular_gcd(y)) return modulargcd(x,y); /* Q[X] */
4034
4035 if (issimplepol(x) || issimplepol(y)) x = RgX_gcd_simple(x,y);
4036 else
4037 {
4038 dx=lg(x); dy=lg(y);
4039 if (dx<dy) { swap(x,y); lswap(dx,dy); }
4040 p1=content(x); p2=content(y); d=ggcd(p1,p2);
4041
4042 tetpil=avma; d = scalarpol(d, vx);
4043 if (dy==3) return gerepile(av,tetpil,d);
4044
4045 av1=avma; lim=stack_lim(av1,1);
4046 u=gdiv(x,p1); v=gdiv(y,p2); g=h=gen_1;
4047 for(;;)
4048 {
4049 GEN r = pseudorem(u,v);
4050 long degq, du, dv, dr=lg(r);
4051
4052 if (dr <= 3)
4053 {
4054 if (gcmp0(r)) break;
4055 avma=av1; return gerepile(av,tetpil,d);
4056 }
4057 if (DEBUGLEVEL > 9) fprintferr("srgcd: dr = %ld\n", dr);
4058 du=lg(u); dv=lg(v); degq=du-dv; u=v;
4059 switch(degq)
4060 {
4061 case 0:
4062 v = gdiv(r,g);
4063 g = leading_term(u);
4064 break;
4065 case 1:
4066 v = gdiv(r,gmul(h,g));
4067 h = g = leading_term(u);
4068 break;
4069 default:
4070 v = gdiv(r,gmul(gpowgs(h,degq),g));
4071 g = leading_term(u);
4072 h = gdiv(gpowgs(g,degq), gpowgs(h,degq-1));
4073 }
4074 if (low_stack(lim, stack_lim(av1,1)))
4075 {
4076 if(DEBUGMEM>1) pari_warn(warnmem,"srgcd");
4077 gerepileall(av1,4,&u,&v,&g,&h);
4078 }
4079 }
4080 p1 = content(v); if (!gcmp1(p1)) v = gdiv(v,p1);
4081 x = gmul(d,v);
4082 }
4083
4084 if (typ(x)!=t_POL) x = scalarpol(x, vx);
4085 else
4086 {
4087 p1=leading_term(x); ty=typ(p1);
4088 if ((ty == t_FRAC || is_intreal_t(ty)) && gsigne(p1)<0) x = gneg(x);
4089 }
4090 return gerepileupto(av,x);
4091 }
4092
4093 GEN
poldisc0(GEN x,long v)4094 poldisc0(GEN x, long v)
4095 {
4096 long tx=typ(x), i;
4097 pari_sp av;
4098 GEN z,p1,p2;
4099
4100 switch(tx)
4101 {
4102 case t_POL:
4103 if (gcmp0(x)) return gen_0;
4104 av = avma; i = 0;
4105 if (v >= 0 && v != varn(x)) x = fix_pol(x,v, &i);
4106 p1 = subres(x, derivpol(x));
4107 p2 = leading_term(x); if (!gcmp1(p2)) p1 = gdiv(p1,p2);
4108 if (degpol(x) & 2) p1 = gneg(p1);
4109 if (i) p1 = gsubst(p1, MAXVARN, pol_x[0]);
4110 return gerepileupto(av,p1);
4111
4112 case t_COMPLEX:
4113 return utoineg(4);
4114
4115 case t_QUAD: case t_POLMOD:
4116 return poldisc0(gel(x,1), v);
4117
4118 case t_QFR: case t_QFI:
4119 av = avma; return gerepileuptoint(av, qf_disc(x));
4120
4121 case t_VEC: case t_COL: case t_MAT:
4122 i=lg(x); z=cgetg(i,tx);
4123 for (i--; i; i--) gel(z,i) = poldisc0(gel(x,i), v);
4124 return z;
4125 }
4126 pari_err(typeer,"discsr");
4127 return NULL; /* not reached */
4128 }
4129
4130 GEN
discsr(GEN x)4131 discsr(GEN x) { return poldisc0(x, -1); }
4132
4133 GEN
reduceddiscsmith(GEN pol)4134 reduceddiscsmith(GEN pol)
4135 {
4136 long i, j, n;
4137 pari_sp av = avma;
4138 GEN polp, p1, m;
4139
4140 if (typ(pol)!=t_POL) pari_err(typeer,"reduceddiscsmith");
4141 n=degpol(pol);
4142 if (n<=0) pari_err(constpoler,"reduceddiscsmith");
4143 check_ZX(pol,"poldiscreduced");
4144 if (!gcmp1(gel(pol,n+2)))
4145 pari_err(talker,"non-monic polynomial in poldiscreduced");
4146 polp = derivpol(pol);
4147 m=cgetg(n+1,t_MAT);
4148 for (j=1; j<=n; j++)
4149 {
4150 p1=cgetg(n+1,t_COL); gel(m,j) = p1;
4151 for (i=1; i<=n; i++) gel(p1,i) = truecoeff(polp,i-1);
4152 if (j<n) polp = grem(RgX_shift_shallow(polp, 1), pol);
4153 }
4154 return gerepileupto(av, smith(m));
4155 }
4156
4157 /***********************************************************************/
4158 /** **/
4159 /** STURM ALGORITHM **/
4160 /** (number of real roots of x in ]a,b]) **/
4161 /** **/
4162 /***********************************************************************/
4163
4164 /* if a (resp. b) is NULL, set it to -oo (resp. +oo) */
4165 long
sturmpart(GEN x,GEN a,GEN b)4166 sturmpart(GEN x, GEN a, GEN b)
4167 {
4168 long sl, sr, s, t, r1;
4169 pari_sp av = avma, lim = stack_lim(av, 1);
4170 GEN g,h,u,v;
4171
4172 if (gcmp0(x)) pari_err(zeropoler,"sturm");
4173 t = typ(x);
4174 if (t != t_POL)
4175 {
4176 if (t == t_INT || t == t_REAL || t == t_FRAC) return 0;
4177 pari_err(typeer,"sturm");
4178 }
4179 s=lg(x); if (s==3) return 0;
4180
4181 sl = gsigne(leading_term(x));
4182 if (s==4)
4183 {
4184 t = a? gsigne(poleval(x,a)): -sl;
4185 if (t == 0) { avma = av; return 0; }
4186 s = b? gsigne(poleval(x,b)): sl;
4187 avma = av; return (s == t)? 0: 1;
4188 }
4189 u=gdiv(x,content(x)); v=derivpol(x);
4190 v=gdiv(v,content(v));
4191 g=gen_1; h=gen_1;
4192 s = b? gsigne(poleval(u,b)): sl;
4193 t = a? gsigne(poleval(u,a)): ((lg(u)&1)? sl: -sl);
4194 r1=0;
4195 sr = b? gsigne(poleval(v,b)): s;
4196 if (sr)
4197 {
4198 if (!s) s=sr;
4199 else if (sr!=s) { s= -s; r1--; }
4200 }
4201 sr = a? gsigne(poleval(v,a)): -t;
4202 if (sr)
4203 {
4204 if (!t) t=sr;
4205 else if (sr!=t) { t= -t; r1++; }
4206 }
4207 for(;;)
4208 {
4209 GEN p1, r = pseudorem(u,v);
4210 long du=lg(u), dv=lg(v), dr=lg(r), degq=du-dv;
4211
4212 if (dr<=2) pari_err(talker,"not a squarefree polynomial in sturm");
4213 if (gsigne(leading_term(v)) > 0 || degq&1) r=gneg_i(r);
4214 sl = gsigne(gel(r,dr-1));
4215 sr = b? gsigne(poleval(r,b)): sl;
4216 if (sr)
4217 {
4218 if (!s) s=sr;
4219 else if (sr!=s) { s= -s; r1--; }
4220 }
4221 sr = a? gsigne(poleval(r,a)): ((dr&1)? sl: -sl);
4222 if (sr)
4223 {
4224 if (!t) t=sr;
4225 else if (sr!=t) { t= -t; r1++; }
4226 }
4227 if (dr==3) { avma=av; return r1; }
4228
4229 u=v; p1 = g; g = gabs(leading_term(u),DEFAULTPREC);
4230 switch(degq)
4231 {
4232 case 0: break;
4233 case 1:
4234 p1 = gmul(h,p1); h = g; break;
4235 default:
4236 p1 = gmul(gpowgs(h,degq),p1);
4237 h = gdivexact(gpowgs(g,degq), gpowgs(h,degq-1));
4238 }
4239 v = gdivexact(r,p1);
4240 if (low_stack(lim,stack_lim(av,1)))
4241 {
4242 if(DEBUGMEM>1) pari_warn(warnmem,"polsturm, dr = %ld",dr);
4243 gerepileall(av,4,&u,&v,&g,&h);
4244 }
4245 }
4246 }
4247
4248 /***********************************************************************/
4249 /** **/
4250 /** GENERIC EXTENDED GCD **/
4251 /** **/
4252 /***********************************************************************/
4253
4254 static GEN
RgXQ_inv_inexact(GEN x,GEN y)4255 RgXQ_inv_inexact(GEN x, GEN y)
4256 {
4257 pari_sp av = avma;
4258 long i, dx = degpol(x), dy = degpol(y), dz = dx+dy;
4259 GEN v, z;
4260
4261 if (dx < 0 || dy < 0) pari_err(talker,"non-invertible polynomial in RgXQ_inv");
4262 v = gauss(sylvestermatrix(y,x), col_ei(dz, dz));
4263 z = cgetg(dy+2,t_POL); z[1] = y[1];
4264 for (i=2; i<dy+2; i++) z[i] = v[dz-i+2];
4265 return gerepilecopy(av, normalizepol_i(z, dy+2));
4266 }
4267 /* assume typ(x) = t_POL */
4268 static GEN
RgXQ_inv(GEN x,GEN y)4269 RgXQ_inv(GEN x, GEN y)
4270 {
4271 long vx=varn(x), vy=varn(y);
4272 pari_sp av, av1;
4273 GEN u, v, d;
4274
4275 while (vx != vy)
4276 {
4277 if (varncmp(vx,vy) > 0)
4278 {
4279 if (vx == BIGINT) return ginv(x);
4280 return gred_rfrac_simple(gen_1, x);
4281 }
4282 if (lg(x)!=3) pari_err(talker,"non-invertible polynomial in RgXQ_inv");
4283 x = gel(x,2); vx = gvar(x);
4284 }
4285 if (isinexact(x) || isinexact(y)) return RgXQ_inv_inexact(x,y);
4286
4287 av = avma; d = subresext(x,y,&u,&v);
4288 if (gcmp0(d)) pari_err(talker,"non-invertible polynomial in RgXQ_inv");
4289 if (typ(d) == t_POL && varn(d) == vx)
4290 {
4291 if (lg(d) > 3) pari_err(talker,"non-invertible polynomial in RgXQ_inv");
4292 d = gel(d,2);
4293 }
4294 av1 = avma; return gerepile(av,av1, gdiv(u,d));
4295 }
4296
4297 GEN
gbezout(GEN x,GEN y,GEN * u,GEN * v)4298 gbezout(GEN x, GEN y, GEN *u, GEN *v)
4299 {
4300 if (typ(x) == t_INT && typ(y) == t_INT) return bezout(x,y,u,v);
4301 return RgX_extgcd(x,y,u,v);
4302 }
4303
4304 GEN
vecbezout(GEN x,GEN y)4305 vecbezout(GEN x, GEN y)
4306 {
4307 GEN z=cgetg(4,t_VEC);
4308 gel(z,3) = gbezout(x,y,(GEN*)(z+1),(GEN*)(z+2));
4309 return z;
4310 }
4311
4312 GEN
vecbezoutres(GEN x,GEN y)4313 vecbezoutres(GEN x, GEN y)
4314 {
4315 GEN z=cgetg(4,t_VEC);
4316 gel(z,3) = subresext(x,y,(GEN*)(z+1),(GEN*)(z+2));
4317 return z;
4318 }
4319
4320
4321 /*******************************************************************/
4322 /* */
4323 /* GENERIC (modular) INVERSE */
4324 /* */
4325 /*******************************************************************/
4326
4327 GEN
ginvmod(GEN x,GEN y)4328 ginvmod(GEN x, GEN y)
4329 {
4330 long tx=typ(x);
4331
4332 switch(typ(y))
4333 {
4334 case t_POL:
4335 if (tx==t_POL) return RgXQ_inv(x,y);
4336 if (is_scalar_t(tx)) return ginv(x);
4337 break;
4338 case t_INT:
4339 if (tx==t_INT) return Fp_inv(x,y);
4340 if (tx==t_POL) return gen_0;
4341 }
4342 pari_err(typeer,"ginvmod");
4343 return NULL; /* not reached */
4344 }
4345
4346 /***********************************************************************/
4347 /** **/
4348 /** NEWTON POLYGON **/
4349 /** **/
4350 /***********************************************************************/
4351
4352 /* assume leading coeff of x is non-zero */
4353 GEN
newtonpoly(GEN x,GEN p)4354 newtonpoly(GEN x, GEN p)
4355 {
4356 GEN y;
4357 long n,ind,a,b,c,u1,u2,r1,r2;
4358 long *vval, num[] = {evaltyp(t_INT) | _evallg(3), 0, 0};
4359
4360 if (typ(x)!=t_POL) pari_err(typeer,"newtonpoly");
4361 n=degpol(x); if (n<=0) { y=cgetg(1,t_VEC); return y; }
4362 y = cgetg(n+1,t_VEC); x += 2; /* now x[i] = term of degree i */
4363 vval = (long *) gpmalloc(sizeof(long)*(n+1));
4364 for (a=0; a<=n; a++) vval[a] = ggval(gel(x,a),p);
4365 for (a=0, ind=1; a<n; a++)
4366 {
4367 if (vval[a] != VERYBIGINT) break;
4368 gel(y,ind++) = utoipos(VERYBIGINT);
4369 }
4370 for (b=a+1; b<=n; a=b, b=a+1)
4371 {
4372 while (vval[b] == VERYBIGINT) b++;
4373 u1=vval[a]-vval[b]; u2=b-a;
4374 for (c=b+1; c<=n; c++)
4375 {
4376 if (vval[c] == VERYBIGINT) continue;
4377 r1=vval[a]-vval[c]; r2=c-a;
4378 if (u1*r2<=u2*r1) { u1=r1; u2=r2; b=c; }
4379 }
4380 while (ind<=b) { affsi(u1,num); gel(y,ind++) = gdivgs(num,u2); }
4381 }
4382 free(vval); return y;
4383 }
4384
4385 static GEN
lift_to_frac(GEN t,GEN mod,GEN amax,GEN bmax,GEN denom)4386 lift_to_frac(GEN t, GEN mod, GEN amax, GEN bmax, GEN denom)
4387 {
4388 GEN a, b;
4389 if (signe(t) < 0) t = addii(t, mod); /* in case t is a centerlift */
4390 if (!ratlift(t, mod, &a, &b, amax, bmax)
4391 || (denom && !dvdii(denom,b))
4392 || !gcmp1(gcdii(a,b))) return NULL;
4393 if (!is_pm1(b)) a = mkfrac(a, b);
4394 return a;
4395 }
4396
4397 /* compute rational lifting for all the components of M modulo mod.
4398 * If one components fails, return NULL.
4399 * See ratlift.
4400 * If denom is not NULL, check that the denominators divide denom
4401 *
4402 * FIXME: NOT stack clean ! a & b stay on the stack.
4403 * If we suppose mod and denom coprime, then a and b are coprime
4404 * so we can do a cgetg(t_FRAC). */
4405 GEN
matratlift(GEN M,GEN mod,GEN amax,GEN bmax,GEN denom)4406 matratlift(GEN M, GEN mod, GEN amax, GEN bmax, GEN denom)
4407 {
4408 pari_sp ltop = avma;
4409 GEN N, a;
4410 long l2, l3;
4411 long i, j;
4412 if (typ(M)!=t_MAT) pari_err(typeer,"matratlift");
4413 l2 = lg(M); l3 = lg(gel(M,1));
4414 N = cgetg(l2, t_MAT);
4415 for (j = 1; j < l2; ++j)
4416 {
4417 gel(N,j) = cgetg(l3, t_COL);
4418 for (i = 1; i < l3; ++i)
4419 {
4420 a = lift_to_frac(gcoeff(M,i,j), mod, amax,bmax,denom);
4421 if (!a) { avma = ltop; return NULL; }
4422 gcoeff(N, i, j) = a;
4423 }
4424 }
4425 return N;
4426 }
4427
4428 GEN
polratlift(GEN P,GEN mod,GEN amax,GEN bmax,GEN denom)4429 polratlift(GEN P, GEN mod, GEN amax, GEN bmax, GEN denom)
4430 {
4431 pari_sp ltop = avma;
4432 GEN Q,a;
4433 long l2, j;
4434 if (typ(P)!=t_POL) pari_err(typeer,"polratlift");
4435 l2 = lg(P);
4436 Q = cgetg(l2, t_POL); Q[1] = P[1];
4437 for (j = 2; j < l2; ++j)
4438 {
4439 a = lift_to_frac(gel(P,j), mod, amax,bmax,denom);
4440 if (!a) { avma = ltop; return NULL; }
4441 gel(Q,j) = a;
4442 }
4443 return Q;
4444 }
4445
4446 /* P,Q in Z[X,Y], nf in Z[Y] irreducible. compute GCD in Q[Y]/(nf)[X].
4447 *
4448 * We essentially follow M. Encarnacion "On a modular Algorithm for computing
4449 * GCDs of polynomials over number fields" (ISSAC'94).
4450 *
4451 * We procede as follows
4452 * 1:compute the gcd modulo primes discarding bad primes as they are detected.
4453 * 2:reconstruct the result via matratlift, stoping as soon as we get weird
4454 * denominators.
4455 * 3:if matratlift succeeds, try the full division.
4456 * Suppose accuracy is insufficient to get the result right: matratlift will
4457 * rarely succeed, and even if it does the polynomial we get has sensible
4458 * coefficients, so the full division will not be too costly.
4459 *
4460 * FIXME: Handle rational coefficient for P and Q.
4461 * If not NULL, den must a a multiple of the denominator of the gcd,
4462 * for example the discriminant of nf.
4463 *
4464 * NOTE: if nf is not irreducible, nfgcd may loop forever, esp. if gcd | nf */
4465 GEN
nfgcd(GEN P,GEN Q,GEN nf,GEN den)4466 nfgcd(GEN P, GEN Q, GEN nf, GEN den)
4467 {
4468 pari_sp ltop = avma;
4469 GEN sol, mod = NULL;
4470 long x = varn(P);
4471 long y = varn(nf);
4472 long d = degpol(nf);
4473 GEN lP, lQ;
4474 if (!signe(P) || !signe(Q)) return zeropol(x);
4475 /*Compute denominators*/
4476 if (!den) den = ZX_disc(nf);
4477 lP = leading_term(P);
4478 lQ = leading_term(Q);
4479 if ( !((typ(lP)==t_INT && is_pm1(lP)) || (typ(lQ)==t_INT && is_pm1(lQ))) )
4480 den = mulii(den, gcdii(ZX_resultant(lP, nf), ZX_resultant(lQ, nf)));
4481 { /*Modular GCD*/
4482 pari_sp btop = avma, st_lim = stack_lim(btop, 1);
4483 ulong p;
4484 long dM=0, dR;
4485 GEN M, dsol;
4486 GEN R, ax, bo;
4487 byteptr primepointer;
4488 for (p = 27449, primepointer = diffptr + 3000; ; )
4489 {
4490 NEXT_PRIME_VIADIFF_CHECK(p, primepointer);
4491 /*Discard primes dividing disc(T) or lc(PQ) */
4492 if (!smodis(den, p)) continue;
4493 if (DEBUGLEVEL>5) fprintferr("nfgcd: p=%d\n",p);
4494 /*Discard primes when modular gcd does not exist*/
4495 if ((R = FlxqX_safegcd(ZXX_to_FlxX(P,p,y),
4496 ZXX_to_FlxX(Q,p,y),
4497 ZX_to_Flx(nf,p), p)) == NULL) continue;
4498 dR = degpol(R);
4499 if (dR == 0) return scalarpol(gen_1, x);
4500 if (mod && dR > dM) continue; /* p divides Res(P/gcd, Q/gcd). Discard. */
4501
4502 R = RgXX_to_RgM(FlxX_to_ZXX(R), d);
4503 /* previous primes divided Res(P/gcd, Q/gcd)? Discard them. */
4504 if (!mod || dR < dM) { M = R; mod = utoipos(p); dM = dR; continue; }
4505 if (low_stack(st_lim, stack_lim(btop, 1)))
4506 {
4507 if (DEBUGMEM>1) pari_warn(warnmem,"nfgcd");
4508 gerepileall(btop, 2, &M, &mod);
4509 }
4510
4511 ax = mulis(Fp_inv(utoipos(p), mod), p);
4512 M = gadd(R, gmul(ax, gsub(M, R)));
4513 mod = mulis(mod, p);
4514 M = lift(FpM_to_mod(M, mod));
4515 /* I suspect it must be better to take amax > bmax*/
4516 bo = sqrti(shifti(mod, -1));
4517 if ((sol = matratlift(M, mod, bo, bo, den)) == NULL) continue;
4518 sol = RgM_to_RgXX(sol,x,y);
4519 dsol = primpart(sol);
4520 if (gcmp0(pseudorem_i(P, dsol, nf))
4521 && gcmp0(pseudorem_i(Q, dsol, nf))) break;
4522 }
4523 }
4524 return gerepilecopy(ltop, sol);
4525 }
4526