Line data Source code
1 : /* Copyright (C) 2000-2004 The PARI group.
2 :
3 : This file is part of the PARI/GP package.
4 :
5 : PARI/GP is free software; you can redistribute it and/or modify it under the
6 : terms of the GNU General Public License as published by the Free Software
7 : Foundation; either version 2 of the License, or (at your option) any later
8 : version. It is distributed in the hope that it will be useful, but WITHOUT
9 : ANY WARRANTY WHATSOEVER.
10 :
11 : Check the License for details. You should have received a copy of it, along
12 : with the package; see the file 'COPYING'. If not, write to the Free Software
13 : Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. */
14 :
15 : /***********************************************************************/
16 : /** **/
17 : /** ARITHMETIC OPERATIONS ON POLYNOMIALS **/
18 : /** (first part) **/
19 : /** **/
20 : /***********************************************************************/
21 : #include "pari.h"
22 : #include "paripriv.h"
23 : /*******************************************************************/
24 : /* */
25 : /* POLYNOMIAL EUCLIDEAN DIVISION */
26 : /* */
27 : /*******************************************************************/
28 : /* x t_POLMOD, y t_POL in the same variable as x[1], return x % y */
29 : static GEN
30 13790 : polmod_mod(GEN x, GEN y)
31 : {
32 13790 : GEN z, a, T = gel(x,1);
33 13790 : if (RgX_equal(T, y)) return gcopy(x);
34 14 : z = cgetg(3,t_POLMOD); T = RgX_gcd(T,y); a = gel(x,2);
35 14 : gel(z,1) = T;
36 14 : gel(z,2) = (typ(a)==t_POL && varn(a)==varn(T))? RgX_rem(a, T): gcopy(a);
37 14 : return z;
38 : }
39 : /* x,y two "scalars", return 0 with type info */
40 : static GEN
41 1575 : rem_scal_scal(GEN x, GEN y)
42 : {
43 1575 : pari_sp av = avma;
44 1575 : GEN z = gadd(gmul(gen_0,x), gmul(gen_0,y));
45 1575 : if (gequal0(y)) pari_err_INV("grem",y);
46 1575 : return gerepileupto(av, simplify(z));
47 : }
48 : /* x pol, y "scalar", return 0 with type info */
49 : static GEN
50 126 : rem_pol_scal(GEN x, GEN y)
51 : {
52 126 : pari_sp av = avma;
53 126 : if (gequal0(y)) pari_err_INV("grem",y);
54 126 : return gerepileupto(av, simplify(gmul(Rg_get_0(x),y)));
55 : }
56 : /* x "scalar", y pol, return x % y with type info */
57 : static GEN
58 1090891 : rem_scal_pol(GEN x, GEN y)
59 : {
60 1090891 : if (degpol(y))
61 : {
62 1089323 : if (!signe(y)) pari_err_INV("grem",y);
63 1089323 : return gmul(x, Rg_get_1(y));
64 : }
65 1568 : y = gel(y,2); return rem_scal_scal(x,y);
66 : }
67 : GEN
68 273 : poldivrem(GEN x, GEN y, GEN *pr)
69 : {
70 273 : const char *f = "euclidean division";
71 273 : long tx = typ(x), ty = typ(y), vx = gvar(x), vy = gvar(y);
72 : GEN z;
73 :
74 273 : if (!is_extscalar_t(tx) || !is_extscalar_t(ty)) pari_err_TYPE2(f,x,y);
75 273 : if (vx == vy && ((tx==t_POLMOD) ^ (ty==t_POLMOD))) pari_err_TYPE2(f,x,y);
76 259 : if (ty != t_POL || varncmp(vx, vy) < 0) /* y "scalar" */
77 : {
78 70 : if (!pr || pr == ONLY_DIVIDES) return gdiv(x,y);
79 70 : if (tx != t_POL || varncmp(vy, vx) < 0) /* x "scalar" */
80 0 : z = rem_scal_scal(x,y);
81 : else
82 70 : z = rem_pol_scal(x,y);
83 70 : if (pr == ONLY_REM) return z;
84 70 : *pr = z; return gdiv(x,y);
85 : }
86 189 : if (tx != t_POL || varncmp(vx, vy) > 0) /* x "scalar" */
87 : {
88 84 : if (!degpol(y)) /* constant t_POL, treat as scalar */
89 : {
90 7 : y = gel(y,2);
91 7 : if (!pr || pr == ONLY_DIVIDES) gdiv(x,y);
92 7 : z = rem_scal_scal(x,y);
93 7 : if (pr == ONLY_REM) return z;
94 7 : *pr = z; return gdiv(x,y);
95 : }
96 77 : if (!signe(y)) pari_err_INV("poldivrem",y);
97 77 : if (!pr || pr == ONLY_DIVIDES) return gequal0(x)? Rg_get_0(y): NULL;
98 77 : z = gmul(x, Rg_get_1(y));
99 77 : if (pr == ONLY_REM) return z;
100 77 : *pr = z; return Rg_get_0(y);
101 : }
102 105 : return RgX_divrem(x,y,pr);
103 : }
104 : GEN
105 637 : gdeuc(GEN x, GEN y)
106 : {
107 637 : const char *f = "euclidean division";
108 637 : long tx = typ(x), ty = typ(y), vx = gvar(x), vy = gvar(y);
109 637 : if (!is_extscalar_t(tx) || !is_extscalar_t(ty)) pari_err_TYPE2(f,x,y);
110 623 : if (vx == vy && ((tx==t_POLMOD) ^ (ty==t_POLMOD))) pari_err_TYPE2(f,x,y);
111 595 : if (ty != t_POL || varncmp(vx, vy) < 0) return gdiv(x,y); /* y "scalar" */
112 455 : if (tx != t_POL || varncmp(vx, vy) > 0)
113 : { /* x "scalar" */
114 140 : if (!signe(y)) pari_err_INV("gdeuc",y);
115 140 : if (!degpol(y)) return gdiv(x, gel(y,2)); /* constant */
116 140 : return Rg_get_0(y);
117 : }
118 315 : return RgX_div(x,y);
119 : }
120 : GEN
121 4295741 : grem(GEN x, GEN y)
122 : {
123 4295741 : const char *f = "euclidean division";
124 4295741 : long tx = typ(x), ty = typ(y), vx = gvar(x), vy = gvar(y);
125 :
126 4295745 : if (ty == t_POL)
127 : {
128 4295682 : if (varncmp(vx,vy) >= 0)
129 : {
130 : pari_sp av;
131 : GEN z;
132 4295681 : if (!signe(y)) pari_err_INV("grem",y);
133 4295682 : if (vx != vy) return rem_scal_pol(x,y);
134 3204791 : switch(tx)
135 : {
136 13790 : case t_POLMOD: return polmod_mod(x,y);
137 3178870 : case t_POL: return RgX_rem(x,y);
138 12082 : case t_RFRAC:
139 : {
140 12082 : GEN a = gel(x,1), b = gel(x,2), p, pol;
141 12082 : if (typ(a) == t_POL && RgX_is_ZX(y) && ZX_is_monic(y))
142 : {
143 12047 : long pa, t = RgX_type2(a,b, &p,&pol,&pa);
144 12047 : if (t == t_FRAC || t == t_INT) return QXQ_div(a, b, y);
145 : }
146 35 : av = avma; z = RgXQ_inv(RgX_rem(b, y), y);
147 28 : return gerepileupto(av, grem(gmul(a, z), y));
148 : }
149 49 : case t_SER:
150 49 : if (RgX_is_monomial(y))
151 : {
152 28 : if (lg(x)-2 + valser(x) < degpol(y)) pari_err_OP("%",x,y);
153 21 : av = avma;
154 21 : return gerepileupto(av, gmod(ser2rfrac_i(x), y));
155 : }
156 21 : default: pari_err_TYPE2("%",x,y);
157 : }
158 : }
159 0 : else switch(tx)
160 : {
161 0 : case t_POL:
162 0 : case t_RFRAC: return rem_pol_scal(x,y);
163 0 : default: pari_err_TYPE2("%",x,y);
164 : }
165 : }
166 63 : if (!is_extscalar_t(tx) || !is_extscalar_t(ty)) pari_err_TYPE2(f,x,y);
167 63 : if (vx == vy && ty==t_POLMOD) pari_err_TYPE2(f,x,y);
168 56 : if (tx != t_POL || varncmp(vx,vy) > 0)
169 : { /* x a "scalar" */
170 0 : if (ty != t_POL || varncmp(vx, vy) < 0) return rem_scal_scal(x,y);
171 0 : return rem_scal_pol(x,y);
172 : }
173 56 : if (ty != t_POL || varncmp(vx, vy) < 0) /* y a "scalar" */
174 56 : return rem_pol_scal(x,y);
175 0 : return RgX_rem(x,y);
176 : }
177 :
178 : /*******************************************************************/
179 : /* */
180 : /* CONVERSIONS RELATED TO p-ADICS */
181 : /* */
182 : /*******************************************************************/
183 : /* x t_PADIC, p a prime or NULL (unset). Consistency check */
184 : static void
185 336 : check_padic_p(GEN x, GEN p)
186 : {
187 336 : GEN q = padic_p(x);
188 336 : if (p && !equalii(p, q)) pari_err_MODULUS("Zp_to_Z", p,q);
189 315 : }
190 : /* shallow */
191 : static GEN
192 4711 : Zp_to_Z(GEN x, GEN p) {
193 4711 : switch(typ(x))
194 : {
195 4473 : case t_INT: break;
196 238 : case t_PADIC:
197 238 : check_padic_p(x, p);
198 217 : x = gtrunc(x); break;
199 0 : default: pari_err_TYPE("Zp_to_Z",x);
200 : }
201 4690 : return x;
202 : }
203 : /* shallow */
204 : static GEN
205 791 : ZpX_to_ZX(GEN x, GEN p)
206 5327 : { pari_APPLY_pol_normalized(Zp_to_Z(gel(x,i), p)); }
207 :
208 : static GEN
209 742 : get_padic_content(GEN f, GEN p)
210 : {
211 742 : GEN c = content(f);
212 742 : if (gequal0(c)) /* O(p^n) can occur */
213 : {
214 0 : if (typ(c) != t_PADIC) pari_err_TYPE("QpX_to_ZX",f);
215 0 : check_padic_p(c, p);
216 0 : c = powis(p, valp(c));
217 : }
218 742 : return c;
219 : }
220 : /* make f suitable for [root|factor]padic. Shallow */
221 : static GEN
222 679 : QpX_to_ZX(GEN f, GEN p)
223 : {
224 679 : GEN c = get_padic_content(f, p);
225 679 : f = RgX_Rg_div(f, c);
226 679 : return ZpX_to_ZX(f, p);
227 : }
228 :
229 : /* x in Z return x + O(pr), pr = p^r. Shallow */
230 : static GEN
231 4816 : Z_to_Zp(GEN x, GEN p, GEN pr, long r)
232 : {
233 4816 : long v, sx = signe(x);
234 :
235 4816 : if (!sx) return zeropadic_shallow(p,r);
236 4193 : v = Z_pvalrem(x,p,&x);
237 4193 : if (v) {
238 973 : if (r <= v) return zeropadic_shallow(p,r);
239 854 : r -= v;
240 854 : pr = powiu(p,r);
241 : }
242 4074 : retmkpadic(modii(x,pr), p, pr, v, r);
243 : }
244 : /* shallow */
245 : static GEN
246 56 : ZV_to_ZpV(GEN z, GEN p, long r)
247 : {
248 56 : long i, l = lg(z);
249 56 : GEN Z = cgetg(l, typ(z)), q = powiu(p, r);
250 161 : for (i=1; i<l; i++) gel(Z,i) = Z_to_Zp(gel(z,i),p,q,r);
251 56 : return Z;
252 : }
253 : /* shallow */
254 : static GEN
255 1372 : ZX_to_ZpX(GEN z, GEN p, GEN q, long r)
256 : {
257 1372 : long i, l = lg(z);
258 1372 : GEN Z = cgetg(l, t_POL); Z[1] = z[1];
259 5880 : for (i=2; i<l; i++) gel(Z,i) = Z_to_Zp(gel(z,i),p,q,r);
260 1372 : return Z;
261 : }
262 : /* return (x + O(p^r)) normalized (multiply by a unit such that leading coeff
263 : * is a power of p), x in Z[X] (or Z_p[X]). Shallow */
264 : static GEN
265 1295 : ZX_to_ZpX_normalized(GEN x, GEN p, GEN pr, long r)
266 : {
267 1295 : long i, lx = lg(x);
268 1295 : GEN z, lead = leading_coeff(x);
269 :
270 1295 : if (gequal1(lead)) return ZX_to_ZpX(x, p, pr, r);
271 63 : (void)Z_pvalrem(lead, p, &lead); lead = Fp_inv(lead, pr);
272 63 : z = cgetg(lx,t_POL);
273 266 : for (i=2; i < lx; i++) gel(z,i) = Z_to_Zp(mulii(lead,gel(x,i)),p,pr,r);
274 63 : z[1] = x[1]; return z;
275 : }
276 : static GEN
277 49 : ZXV_to_ZpXQV(GEN z, GEN T, GEN p, long r)
278 : {
279 49 : long i, l = lg(z);
280 49 : GEN Z = cgetg(l, typ(z)), q = powiu(p, r);
281 49 : T = ZX_copy(T);
282 126 : for (i=1; i<lg(z); i++) gel(Z,i) = mkpolmod(ZX_to_ZpX(gel(z,i),p,q,r),T);
283 49 : return Z;
284 : }
285 : /* shallow */
286 : static GEN
287 63 : QpXQX_to_ZXY(GEN f, GEN p)
288 : {
289 63 : GEN c = get_padic_content(f, p);
290 63 : long i, l = lg(f);
291 63 : f = RgX_Rg_div(f,c);
292 287 : for (i=2; i<l; i++)
293 : {
294 231 : GEN t = gel(f,i);
295 231 : switch(typ(t))
296 : {
297 77 : case t_POLMOD:
298 77 : t = gel(t,2);
299 77 : t = (typ(t) == t_POL)? ZpX_to_ZX(t, p): Zp_to_Z(t, p);
300 77 : break;
301 0 : case t_POL: t = ZpX_to_ZX(t, p); break;
302 154 : default: t = Zp_to_Z(t, p); break;
303 : }
304 224 : gel(f,i) = t;
305 : }
306 56 : return f;
307 : }
308 :
309 : /*******************************************************************/
310 : /* */
311 : /* p-ADIC ROOTS */
312 : /* */
313 : /*******************************************************************/
314 :
315 : /* f primitive ZX, squarefree, leading term prime to p; 0 <= a < p such that
316 : * f(a) = 0 mod p. Return p-adic roots of f equal to a mod p, in
317 : * precision >= prec */
318 : GEN
319 2863 : ZX_Zp_root(GEN f, GEN a, GEN p, long prec)
320 : {
321 : GEN z, R;
322 : long i, j, k;
323 :
324 2863 : if (signe(FpX_eval(FpX_deriv(f, p), a, p)))
325 : { /* simple zero mod p, go all the way to p^prec */
326 2632 : if (prec > 1) a = ZpX_liftroot(f, a, p, prec);
327 2632 : return mkcol(a);
328 : }
329 :
330 231 : f = ZX_unscale_div(ZX_translate(f,a), p); /* f(pX + a) / p */
331 231 : (void)ZX_pvalrem(f,p,&f);
332 231 : z = cgetg(degpol(f)+1,t_COL);
333 :
334 231 : R = FpX_roots(f, p);
335 574 : for (j=i=1; i<lg(R); i++)
336 : {
337 343 : GEN u = ZX_Zp_root(f, gel(R,i), p, prec-1);
338 756 : for (k=1; k<lg(u); k++) gel(z,j++) = addii(a, mulii(p, gel(u,k)));
339 : }
340 231 : setlg(z,j); return z;
341 : }
342 :
343 : /* a t_PADIC, return vector of p-adic roots of f equal to a (mod p) */
344 : GEN
345 56 : Zp_appr(GEN f, GEN a)
346 : {
347 56 : pari_sp av = avma;
348 56 : GEN z, p = padic_p(a);
349 56 : long v = valp(a), prec = v;
350 :
351 56 : if (signe(padic_u(a))) prec += precp(a);
352 56 : f = QpX_to_ZX(f, p);
353 42 : if (degpol(f) <= 0) pari_err_CONSTPOL("Zp_appr");
354 42 : if (v < 0) pari_err_DOMAIN("padicappr", "v(a)", "<", gen_0, stoi(v));
355 35 : f = ZX_radical(f);
356 35 : a = padic_to_Fp(a, p);
357 35 : if (signe(FpX_eval(f,a,p))) { set_avma(av); return cgetg(1,t_COL); }
358 28 : z = ZX_Zp_root(f, a, p, prec);
359 28 : return gerepilecopy(av, ZV_to_ZpV(z, p, prec));
360 : }
361 : static long
362 126 : pval(GEN x, GEN p) { return typ(x) == t_INT? Z_pval(x,p): ZX_pval(x,p); }
363 : /* f a ZXX, f(0) != 0 */
364 : static GEN
365 595 : pnormalize(GEN f, GEN p, GEN T, long prec, long n,
366 : GEN *plead, long *pprec, int *prev)
367 : {
368 595 : *plead = leading_coeff(f);
369 595 : *pprec = prec;
370 595 : *prev = 0;
371 595 : if (!isint1(*plead))
372 : {
373 63 : long v = pval(*plead,p), v1 = pval(constant_coeff(f),p);
374 63 : if (v1 < v)
375 : {
376 49 : *prev = 1;
377 49 : f = RgX_recip_i(f); /* f(0) != 0 so degree is the same */
378 : /* beware loss of precision from lc(factor), whose valuation is <= v */
379 49 : *pprec += v; v = v1;
380 : }
381 63 : *pprec += v * n;
382 : }
383 595 : if (!T) return ZX_Q_normalize(f, plead);
384 14 : *plead = gen_1;
385 14 : return FpXQX_normalize(f, T, powiu(p,*pprec));
386 : }
387 :
388 : /**************************************************************************/
389 :
390 : static void
391 238 : scalar_getprec(GEN x, long *pprec, GEN *pp)
392 : {
393 238 : if (typ(x)==t_PADIC)
394 : {
395 98 : long e = valp(x); if (signe(padic_u(x))) e += precp(x);
396 98 : if (e < *pprec) *pprec = e;
397 98 : check_padic_p(x, *pp);
398 98 : *pp = padic_p(x);
399 : }
400 238 : }
401 : static void
402 98 : getprec(GEN x, long *pprec, GEN *pp)
403 : {
404 : long i;
405 98 : if (typ(x) != t_POL) scalar_getprec(x, pprec, pp);
406 : else
407 266 : for (i = lg(x)-1; i>1; i--) scalar_getprec(gel(x,i), pprec, pp);
408 98 : }
409 :
410 : /* assume f(a) = 0 (mod T,p) */
411 : static GEN
412 105 : ZXY_ZpQ_root(GEN f, GEN a, GEN T, GEN p, long prec)
413 : {
414 : GEN z, R;
415 : long i, j, k, lR;
416 105 : if (signe(FqX_eval(FqX_deriv(f,T,p), a, T,p)))
417 : { /* simple zero mod (T,p), go all the way to p^prec */
418 77 : if (prec > 1) a = ZpXQX_liftroot(f, a, T, p, prec);
419 77 : return mkcol(a);
420 : }
421 28 : f = RgX_unscale(RgXQX_translate(f, a, T), p);
422 28 : f = RgX_Rg_div(f, powiu(p, gvaluation(f,p)));
423 28 : z = cgetg(degpol(f)+1,t_COL);
424 28 : R = FpXQX_roots(FqX_red(f,T,p), T, p); lR = lg(R);
425 70 : for(j=i=1; i<lR; i++)
426 : {
427 42 : GEN u = ZXY_ZpQ_root(f, gel(R,i), T, p, prec-1);
428 84 : for (k=1; k<lg(u); k++) gel(z,j++) = gadd(a, gmul(p, gel(u,k)));
429 : }
430 28 : setlg(z,j); return z;
431 : }
432 :
433 : /* a belongs to finite extension of Q_p, return all roots of f equal to a
434 : * mod p. Don't assume f(a) = 0 (mod p) */
435 : GEN
436 105 : padicappr(GEN f, GEN a)
437 : {
438 : GEN p, z, T, Tp;
439 : long prec;
440 105 : pari_sp av = avma;
441 :
442 105 : if (typ(f)!=t_POL) pari_err_TYPE("padicappr",f);
443 105 : switch(typ(a)) {
444 56 : case t_PADIC: return Zp_appr(f,a);
445 49 : case t_POLMOD: break;
446 0 : default: pari_err_TYPE("padicappr",a);
447 : }
448 49 : if (gequal0(f)) pari_err_ROOTS0("padicappr");
449 49 : T = gel(a,1);
450 49 : a = gel(a,2);
451 49 : p = NULL; prec = LONG_MAX;
452 49 : getprec(a, &prec, &p);
453 49 : getprec(T, &prec, &p); if (!p) pari_err_TYPE("padicappr",T);
454 49 : f = QpXQX_to_ZXY(f, p);
455 42 : if (typ(a) != t_POL) a = scalarpol_shallow(a, varn(T));
456 42 : a = ZpX_to_ZX(a,p);
457 42 : T = QpX_to_ZX(T,p);
458 : /* ensure that f /= (f,f') is separable */
459 42 : (void)nfgcd_all(f, RgX_deriv(f), T, NULL, &f);
460 :
461 42 : Tp = FpX_red(T, p); a = FqX_red(a, Tp, p);
462 42 : if (!gequal0(FqX_eval(FqX_red(f,Tp,p), a, Tp,p)))
463 7 : { set_avma(av); return cgetg(1,t_COL); } /* f(a) != 0 (mod p,T) */
464 35 : z = ZXY_ZpQ_root(f, a, T, p, prec);
465 35 : return gerepilecopy(av, ZXV_to_ZpXQV(z, T, p, prec));
466 : }
467 :
468 : /* vector of p-adic roots of the ZX f, leading term prime to p. Shallow */
469 : static GEN
470 35 : ZX_Zp_roots(GEN f, GEN p, long prec)
471 : {
472 : long l, i;
473 : GEN r;
474 :
475 35 : f = ZX_radical(f);
476 35 : r = FpX_roots(f, p);
477 35 : l = lg(r); if (l == 1) return r;
478 91 : for (i = 1; i < l; i++) gel(r,i) = ZX_Zp_root(f, gel(r,i), p, prec);
479 28 : return ZV_to_ZpV(shallowconcat1(r), p, prec);
480 : }
481 : /* vector of p-adic roots of the ZXX f, leading term prime to p. Shallow */
482 : static GEN
483 14 : ZXY_ZpQ_roots(GEN f, GEN T, GEN p, long prec)
484 : {
485 : long l, i;
486 : GEN r;
487 :
488 14 : (void)nfgcd_all(f, RgX_deriv(f), T, NULL, &f);
489 14 : r = FqX_roots(f, FpX_red(T,p), p);
490 14 : l = lg(r); if (l == 1) return r;
491 42 : for (i = 1; i < l; i++) gel(r,i) = ZXY_ZpQ_root(f, gel(r,i), T, p, prec);
492 14 : return ZXV_to_ZpXQV(shallowconcat1(r), T, p, prec);
493 : }
494 :
495 : /* return p-adic roots of f, precision prec */
496 : GEN
497 56 : polrootspadic(GEN f, GEN Tp, long prec)
498 : {
499 56 : pari_sp av = avma;
500 : GEN lead, y, T, p;
501 : long PREC, i, k, v;
502 : int reverse;
503 :
504 56 : if (!ff_parse_Tp(Tp, &T,&p,0)) pari_err_TYPE("polrootspadic",Tp);
505 56 : if (typ(f)!=t_POL) pari_err_TYPE("polrootspadic",f);
506 56 : if (gequal0(f)) pari_err_ROOTS0("polrootspadic");
507 56 : if (prec <= 0)
508 7 : pari_err_DOMAIN("polrootspadic", "precision", "<=",gen_0,stoi(prec));
509 49 : f = T? QpXQX_to_ZXY(f, p): QpX_to_ZX(f, p);
510 49 : v = RgX_valrem(f, &f);
511 49 : f = pnormalize(f, p, T, prec, 1, &lead, &PREC, &reverse);
512 49 : y = T? ZXY_ZpQ_roots(f,T,p,PREC): ZX_Zp_roots(f,p,PREC);
513 49 : k = lg(y);
514 49 : if (lead != gen_1) y = RgC_Rg_div(y, lead);
515 49 : if (reverse)
516 7 : for (i=1; i<k; i++) gel(y,i) = ginv(gel(y,i));
517 49 : if (v) y = shallowconcat(zeropadic_shallow(p, prec), y);
518 49 : return gerepilecopy(av, y);
519 : }
520 :
521 : /*******************************************************************/
522 : /* */
523 : /* FACTORIZATION in Zp[X], using ROUND4 */
524 : /* */
525 : /*******************************************************************/
526 :
527 : int
528 3070 : cmp_padic(GEN x, GEN y)
529 : {
530 : long vx, vy;
531 3070 : if (x == gen_0) return -1;
532 3070 : if (y == gen_0) return 1;
533 3070 : vx = valp(x);
534 3070 : vy = valp(y);
535 3070 : if (vx < vy) return 1;
536 3035 : if (vx > vy) return -1;
537 2776 : return cmpii(padic_u(x), padic_u(y));
538 : }
539 :
540 : /* replace p^e by p*...*p [ factors are not known to be equal, only close at
541 : * input accuracy ] */
542 : static GEN
543 49 : famat_flatten(GEN fa)
544 : {
545 49 : GEN y, P = gel(fa,1), E = gel(fa,2);
546 49 : long i, l = lg(E);
547 49 : y = cgetg(l, t_VEC);
548 147 : for (i = 1; i < l; i++)
549 : {
550 98 : GEN p = gel(P,i);
551 98 : long e = itou(gel(E,i));
552 98 : gel(y,i) = const_col(e, p);
553 : }
554 49 : y = shallowconcat1(y); return mkmat2(y, const_col(lg(y)-1, gen_1));
555 : }
556 :
557 : GEN
558 581 : factorpadic(GEN f, GEN p, long r)
559 : {
560 581 : pari_sp av = avma;
561 : GEN y, ppow;
562 : long v, n;
563 581 : int reverse = 0, exact;
564 :
565 581 : if (typ(f)!=t_POL) pari_err_TYPE("factorpadic",f);
566 581 : if (typ(p)!=t_INT) pari_err_TYPE("factorpadic",p);
567 581 : if (r <= 0) pari_err_DOMAIN("factorpadic", "precision", "<=",gen_0,stoi(r));
568 574 : if (!signe(f)) return prime_fact(f);
569 574 : if (!degpol(f)) return trivial_fact();
570 :
571 574 : exact = RgX_is_QX(f); /* before RgX_valrem which may lose type information */
572 574 : v = RgX_valrem_inexact(f, &f);
573 574 : ppow = powiu(p,r);
574 574 : n = degpol(f);
575 574 : if (!n)
576 28 : y = trivial_fact();
577 : else
578 : {
579 : GEN P, lead;
580 : long i, l, pr;
581 :
582 546 : f = QpX_to_ZX(f, p);
583 546 : f = pnormalize(f, p, NULL, r, n-1, &lead, &pr, &reverse);
584 546 : y = ZpX_monic_factor(f, p, pr);
585 546 : P = gel(y,1); l = lg(P);
586 546 : if (!isint1(lead))
587 294 : for (i=1; i<l; i++) gel(P,i) = Q_primpart(RgX_unscale(gel(P,i), lead));
588 1841 : for (i=1; i<l; i++)
589 : {
590 1295 : GEN t = gel(P,i);
591 1295 : if (reverse) t = RgX_recip_shallow(t);
592 1295 : gel(P,i) = ZX_to_ZpX_normalized(t,p,ppow,r);
593 : }
594 : }
595 574 : if (v)
596 : { /* v > 0 */
597 63 : GEN X = ZX_to_ZpX(pol_x(varn(f)), p, ppow, r);
598 63 : y = famat_mulpow_shallow(y, X, utoipos(v));
599 : }
600 574 : if (!exact) y = famat_flatten(y);
601 574 : return gerepilecopy(av, sort_factor_pol(y, cmp_padic));
602 : }
|
