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 11820 : polmod_mod(GEN x, GEN y)
31 : {
32 11820 : GEN z, a, T = gel(x,1);
33 11820 : if (RgX_equal(T, y)) return gcopy(x);
34 12 : z = cgetg(3,t_POLMOD); T = RgX_gcd(T,y); a = gel(x,2);
35 12 : gel(z,1) = T;
36 12 : gel(z,2) = (typ(a)==t_POL && varn(a)==varn(T))? RgX_rem(a, T): gcopy(a);
37 12 : return z;
38 : }
39 : /* x,y two "scalars", return 0 with type info */
40 : static GEN
41 1350 : rem_scal_scal(GEN x, GEN y)
42 : {
43 1350 : pari_sp av = avma;
44 1350 : GEN z = gadd(gmul(gen_0,x), gmul(gen_0,y));
45 1350 : if (gequal0(y)) pari_err_INV("grem",y);
46 1350 : return gerepileupto(av, simplify(z));
47 : }
48 : /* x pol, y "scalar", return 0 with type info */
49 : static GEN
50 108 : rem_pol_scal(GEN x, GEN y)
51 : {
52 108 : pari_sp av = avma;
53 108 : if (gequal0(y)) pari_err_INV("grem",y);
54 108 : 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 935087 : rem_scal_pol(GEN x, GEN y)
59 : {
60 935087 : if (degpol(y))
61 : {
62 933743 : if (!signe(y)) pari_err_INV("grem",y);
63 933743 : return gmul(x, Rg_get_1(y));
64 : }
65 1344 : y = gel(y,2); return rem_scal_scal(x,y);
66 : }
67 : GEN
68 234 : poldivrem(GEN x, GEN y, GEN *pr)
69 : {
70 234 : const char *f = "euclidean division";
71 234 : long tx = typ(x), ty = typ(y), vx = gvar(x), vy = gvar(y);
72 : GEN z;
73 :
74 234 : if (!is_extscalar_t(tx) || !is_extscalar_t(ty)) pari_err_TYPE2(f,x,y);
75 234 : if (vx == vy && ((tx==t_POLMOD) ^ (ty==t_POLMOD))) pari_err_TYPE2(f,x,y);
76 222 : if (ty != t_POL || varncmp(vx, vy) < 0) /* y "scalar" */
77 : {
78 60 : if (!pr || pr == ONLY_DIVIDES) return gdiv(x,y);
79 60 : if (tx != t_POL || varncmp(vy, vx) < 0) /* x "scalar" */
80 0 : z = rem_scal_scal(x,y);
81 : else
82 60 : z = rem_pol_scal(x,y);
83 60 : if (pr == ONLY_REM) return z;
84 60 : *pr = z; return gdiv(x,y);
85 : }
86 162 : if (tx != t_POL || varncmp(vx, vy) > 0) /* x "scalar" */
87 : {
88 72 : if (!degpol(y)) /* constant t_POL, treat as scalar */
89 : {
90 6 : y = gel(y,2);
91 6 : if (!pr || pr == ONLY_DIVIDES) gdiv(x,y);
92 6 : z = rem_scal_scal(x,y);
93 6 : if (pr == ONLY_REM) return z;
94 6 : *pr = z; return gdiv(x,y);
95 : }
96 66 : if (!signe(y)) pari_err_INV("poldivrem",y);
97 66 : if (!pr || pr == ONLY_DIVIDES) return gequal0(x)? Rg_get_0(y): NULL;
98 66 : z = gmul(x, Rg_get_1(y));
99 66 : if (pr == ONLY_REM) return z;
100 66 : *pr = z; return Rg_get_0(y);
101 : }
102 90 : return RgX_divrem(x,y,pr);
103 : }
104 : GEN
105 546 : gdeuc(GEN x, GEN y)
106 : {
107 546 : const char *f = "euclidean division";
108 546 : long tx = typ(x), ty = typ(y), vx = gvar(x), vy = gvar(y);
109 546 : if (!is_extscalar_t(tx) || !is_extscalar_t(ty)) pari_err_TYPE2(f,x,y);
110 534 : if (vx == vy && ((tx==t_POLMOD) ^ (ty==t_POLMOD))) pari_err_TYPE2(f,x,y);
111 510 : if (ty != t_POL || varncmp(vx, vy) < 0) return gdiv(x,y); /* y "scalar" */
112 390 : if (tx != t_POL || varncmp(vx, vy) > 0)
113 : { /* x "scalar" */
114 120 : if (!signe(y)) pari_err_INV("gdeuc",y);
115 120 : if (!degpol(y)) return gdiv(x, gel(y,2)); /* constant */
116 120 : return Rg_get_0(y);
117 : }
118 270 : return RgX_div(x,y);
119 : }
120 : GEN
121 3681774 : grem(GEN x, GEN y)
122 : {
123 3681774 : const char *f = "euclidean division";
124 3681774 : long tx = typ(x), ty = typ(y), vx = gvar(x), vy = gvar(y);
125 :
126 3681774 : if (ty == t_POL)
127 : {
128 3681720 : if (varncmp(vx,vy) >= 0)
129 : {
130 : pari_sp av;
131 : GEN z;
132 3681720 : if (!signe(y)) pari_err_INV("grem",y);
133 3681720 : if (vx != vy) return rem_scal_pol(x,y);
134 2746633 : switch(tx)
135 : {
136 11820 : case t_POLMOD: return polmod_mod(x,y);
137 2724415 : case t_POL: return RgX_rem(x,y);
138 10356 : case t_RFRAC:
139 : {
140 10356 : GEN a = gel(x,1), b = gel(x,2), p, pol;
141 10356 : if (typ(a) == t_POL && RgX_is_ZX(y) && ZX_is_monic(y))
142 : {
143 10326 : long pa, t = RgX_type2(a,b, &p,&pol,&pa);
144 10326 : if (t == t_FRAC || t == t_INT) return QXQ_div(a, b, y);
145 : }
146 30 : av = avma; z = RgXQ_inv(RgX_rem(b, y), y);
147 24 : return gerepileupto(av, grem(gmul(a, z), y));
148 : }
149 42 : case t_SER:
150 42 : if (RgX_is_monomial(y))
151 : {
152 24 : if (lg(x)-2 + valser(x) < degpol(y)) pari_err_OP("%",x,y);
153 18 : av = avma;
154 18 : return gerepileupto(av, gmod(ser2rfrac_i(x), y));
155 : }
156 18 : 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 54 : if (!is_extscalar_t(tx) || !is_extscalar_t(ty)) pari_err_TYPE2(f,x,y);
167 54 : if (vx == vy && ty==t_POLMOD) pari_err_TYPE2(f,x,y);
168 48 : 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 48 : if (ty != t_POL || varncmp(vx, vy) < 0) /* y a "scalar" */
174 48 : 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 288 : check_padic_p(GEN x, GEN p)
186 : {
187 288 : GEN q = padic_p(x);
188 288 : if (p && !equalii(p, q)) pari_err_MODULUS("Zp_to_Z", p,q);
189 270 : }
190 : /* shallow */
191 : static GEN
192 4038 : Zp_to_Z(GEN x, GEN p) {
193 4038 : switch(typ(x))
194 : {
195 3834 : case t_INT: break;
196 204 : case t_PADIC:
197 204 : check_padic_p(x, p);
198 186 : x = gtrunc(x); break;
199 0 : default: pari_err_TYPE("Zp_to_Z",x);
200 : }
201 4020 : return x;
202 : }
203 : /* shallow */
204 : static GEN
205 678 : ZpX_to_ZX(GEN x, GEN p)
206 4566 : { pari_APPLY_pol_normalized(Zp_to_Z(gel(x,i), p)); }
207 :
208 : static GEN
209 636 : get_padic_content(GEN f, GEN p)
210 : {
211 636 : GEN c = content(f);
212 636 : 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 636 : return c;
219 : }
220 : /* make f suitable for [root|factor]padic. Shallow */
221 : static GEN
222 582 : QpX_to_ZX(GEN f, GEN p)
223 : {
224 582 : GEN c = get_padic_content(f, p);
225 582 : f = RgX_Rg_div(f, c);
226 582 : return ZpX_to_ZX(f, p);
227 : }
228 :
229 : /* x in Z return x + O(pr), pr = p^r. Shallow */
230 : static GEN
231 4128 : Z_to_Zp(GEN x, GEN p, GEN pr, long r)
232 : {
233 4128 : long v, sx = signe(x);
234 :
235 4128 : if (!sx) return zeropadic_shallow(p,r);
236 3594 : v = Z_pvalrem(x,p,&x);
237 3594 : if (v) {
238 834 : if (r <= v) return zeropadic_shallow(p,r);
239 732 : r -= v;
240 732 : pr = powiu(p,r);
241 : }
242 3492 : retmkpadic(modii(x,pr), p, pr, v, r);
243 : }
244 : /* shallow */
245 : static GEN
246 48 : ZV_to_ZpV(GEN z, GEN p, long r)
247 : {
248 48 : long i, l = lg(z);
249 48 : GEN Z = cgetg(l, typ(z)), q = powiu(p, r);
250 138 : for (i=1; i<l; i++) gel(Z,i) = Z_to_Zp(gel(z,i),p,q,r);
251 48 : return Z;
252 : }
253 : /* shallow */
254 : static GEN
255 1176 : ZX_to_ZpX(GEN z, GEN p, GEN q, long r)
256 : {
257 1176 : long i, l = lg(z);
258 1176 : GEN Z = cgetg(l, t_POL); Z[1] = z[1];
259 5040 : for (i=2; i<l; i++) gel(Z,i) = Z_to_Zp(gel(z,i),p,q,r);
260 1176 : 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 1110 : ZX_to_ZpX_normalized(GEN x, GEN p, GEN pr, long r)
266 : {
267 1110 : long i, lx = lg(x);
268 1110 : GEN z, lead = leading_coeff(x);
269 :
270 1110 : if (gequal1(lead)) return ZX_to_ZpX(x, p, pr, r);
271 54 : (void)Z_pvalrem(lead, p, &lead); lead = Fp_inv(lead, pr);
272 54 : z = cgetg(lx,t_POL);
273 228 : for (i=2; i < lx; i++) gel(z,i) = Z_to_Zp(mulii(lead,gel(x,i)),p,pr,r);
274 54 : z[1] = x[1]; return z;
275 : }
276 : static GEN
277 42 : ZXV_to_ZpXQV(GEN z, GEN T, GEN p, long r)
278 : {
279 42 : long i, l = lg(z);
280 42 : GEN Z = cgetg(l, typ(z)), q = powiu(p, r);
281 42 : T = ZX_copy(T);
282 108 : for (i=1; i<lg(z); i++) gel(Z,i) = mkpolmod(ZX_to_ZpX(gel(z,i),p,q,r),T);
283 42 : return Z;
284 : }
285 : /* shallow */
286 : static GEN
287 54 : QpXQX_to_ZXY(GEN f, GEN p)
288 : {
289 54 : GEN c = get_padic_content(f, p);
290 54 : long i, l = lg(f);
291 54 : f = RgX_Rg_div(f,c);
292 246 : for (i=2; i<l; i++)
293 : {
294 198 : GEN t = gel(f,i);
295 198 : switch(typ(t))
296 : {
297 66 : case t_POLMOD:
298 66 : t = gel(t,2);
299 66 : t = (typ(t) == t_POL)? ZpX_to_ZX(t, p): Zp_to_Z(t, p);
300 66 : break;
301 0 : case t_POL: t = ZpX_to_ZX(t, p); break;
302 132 : default: t = Zp_to_Z(t, p); break;
303 : }
304 192 : gel(f,i) = t;
305 : }
306 48 : 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 2454 : ZX_Zp_root(GEN f, GEN a, GEN p, long prec)
320 : {
321 : GEN z, R;
322 : long i, j, k;
323 :
324 2454 : if (signe(FpX_eval(FpX_deriv(f, p), a, p)))
325 : { /* simple zero mod p, go all the way to p^prec */
326 2256 : if (prec > 1) a = ZpX_liftroot(f, a, p, prec);
327 2256 : return mkcol(a);
328 : }
329 :
330 198 : f = ZX_unscale_div(ZX_translate(f,a), p); /* f(pX + a) / p */
331 198 : (void)ZX_pvalrem(f,p,&f);
332 198 : z = cgetg(degpol(f)+1,t_COL);
333 :
334 198 : R = FpX_roots(f, p);
335 492 : for (j=i=1; i<lg(R); i++)
336 : {
337 294 : GEN u = ZX_Zp_root(f, gel(R,i), p, prec-1);
338 648 : for (k=1; k<lg(u); k++) gel(z,j++) = addii(a, mulii(p, gel(u,k)));
339 : }
340 198 : 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 48 : Zp_appr(GEN f, GEN a)
346 : {
347 48 : pari_sp av = avma;
348 48 : GEN z, p = padic_p(a);
349 48 : long v = valp(a), prec = v;
350 :
351 48 : if (signe(padic_u(a))) prec += precp(a);
352 48 : f = QpX_to_ZX(f, p);
353 36 : if (degpol(f) <= 0) pari_err_CONSTPOL("Zp_appr");
354 36 : if (v < 0) pari_err_DOMAIN("padicappr", "v(a)", "<", gen_0, stoi(v));
355 30 : f = ZX_radical(f);
356 30 : a = padic_to_Fp(a, p);
357 30 : if (signe(FpX_eval(f,a,p))) { set_avma(av); return cgetg(1,t_COL); }
358 24 : z = ZX_Zp_root(f, a, p, prec);
359 24 : return gerepilecopy(av, ZV_to_ZpV(z, p, prec));
360 : }
361 : static long
362 108 : 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 510 : pnormalize(GEN f, GEN p, GEN T, long prec, long n,
366 : GEN *plead, long *pprec, int *prev)
367 : {
368 510 : *plead = leading_coeff(f);
369 510 : *pprec = prec;
370 510 : *prev = 0;
371 510 : if (!isint1(*plead))
372 : {
373 54 : long v = pval(*plead,p), v1 = pval(constant_coeff(f),p);
374 54 : if (v1 < v)
375 : {
376 42 : *prev = 1;
377 42 : 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 42 : *pprec += v; v = v1;
380 : }
381 54 : *pprec += v * n;
382 : }
383 510 : if (!T) return ZX_Q_normalize(f, plead);
384 12 : *plead = gen_1;
385 12 : return FpXQX_normalize(f, T, powiu(p,*pprec));
386 : }
387 :
388 : /**************************************************************************/
389 :
390 : static void
391 204 : scalar_getprec(GEN x, long *pprec, GEN *pp)
392 : {
393 204 : if (typ(x)==t_PADIC)
394 : {
395 84 : long e = valp(x); if (signe(padic_u(x))) e += precp(x);
396 84 : if (e < *pprec) *pprec = e;
397 84 : check_padic_p(x, *pp);
398 84 : *pp = padic_p(x);
399 : }
400 204 : }
401 : static void
402 84 : getprec(GEN x, long *pprec, GEN *pp)
403 : {
404 : long i;
405 84 : if (typ(x) != t_POL) scalar_getprec(x, pprec, pp);
406 : else
407 228 : for (i = lg(x)-1; i>1; i--) scalar_getprec(gel(x,i), pprec, pp);
408 84 : }
409 :
410 : /* assume f(a) = 0 (mod T,p) */
411 : static GEN
412 90 : 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 90 : 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 66 : if (prec > 1) a = ZpXQX_liftroot(f, a, T, p, prec);
419 66 : return mkcol(a);
420 : }
421 24 : f = RgX_unscale(RgXQX_translate(f, a, T), p);
422 24 : f = RgX_Rg_div(f, powiu(p, gvaluation(f,p)));
423 24 : z = cgetg(degpol(f)+1,t_COL);
424 24 : R = FpXQX_roots(FqX_red(f,T,p), T, p); lR = lg(R);
425 60 : for(j=i=1; i<lR; i++)
426 : {
427 36 : GEN u = ZXY_ZpQ_root(f, gel(R,i), T, p, prec-1);
428 72 : for (k=1; k<lg(u); k++) gel(z,j++) = gadd(a, gmul(p, gel(u,k)));
429 : }
430 24 : 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 90 : padicappr(GEN f, GEN a)
437 : {
438 : GEN p, z, T, Tp;
439 : long prec;
440 90 : pari_sp av = avma;
441 :
442 90 : if (typ(f)!=t_POL) pari_err_TYPE("padicappr",f);
443 90 : switch(typ(a)) {
444 48 : case t_PADIC: return Zp_appr(f,a);
445 42 : case t_POLMOD: break;
446 0 : default: pari_err_TYPE("padicappr",a);
447 : }
448 42 : if (gequal0(f)) pari_err_ROOTS0("padicappr");
449 42 : T = gel(a,1);
450 42 : a = gel(a,2);
451 42 : p = NULL; prec = LONG_MAX;
452 42 : getprec(a, &prec, &p);
453 42 : getprec(T, &prec, &p); if (!p) pari_err_TYPE("padicappr",T);
454 42 : f = QpXQX_to_ZXY(f, p);
455 36 : if (typ(a) != t_POL) a = scalarpol_shallow(a, varn(T));
456 36 : a = ZpX_to_ZX(a,p);
457 36 : T = QpX_to_ZX(T,p);
458 : /* ensure that f /= (f,f') is separable */
459 36 : (void)nfgcd_all(f, RgX_deriv(f), T, NULL, &f);
460 :
461 36 : Tp = FpX_red(T, p); a = FqX_red(a, Tp, p);
462 36 : if (!gequal0(FqX_eval(FqX_red(f,Tp,p), a, Tp,p)))
463 6 : { set_avma(av); return cgetg(1,t_COL); } /* f(a) != 0 (mod p,T) */
464 30 : z = ZXY_ZpQ_root(f, a, T, p, prec);
465 30 : 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 30 : ZX_Zp_roots(GEN f, GEN p, long prec)
471 : {
472 : long l, i;
473 : GEN r;
474 :
475 30 : f = ZX_radical(f);
476 30 : r = FpX_roots(f, p);
477 30 : l = lg(r); if (l == 1) return r;
478 78 : for (i = 1; i < l; i++) gel(r,i) = ZX_Zp_root(f, gel(r,i), p, prec);
479 24 : 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 12 : ZXY_ZpQ_roots(GEN f, GEN T, GEN p, long prec)
484 : {
485 : long l, i;
486 : GEN r;
487 :
488 12 : (void)nfgcd_all(f, RgX_deriv(f), T, NULL, &f);
489 12 : r = FqX_roots(f, FpX_red(T,p), p);
490 12 : l = lg(r); if (l == 1) return r;
491 36 : for (i = 1; i < l; i++) gel(r,i) = ZXY_ZpQ_root(f, gel(r,i), T, p, prec);
492 12 : return ZXV_to_ZpXQV(shallowconcat1(r), T, p, prec);
493 : }
494 :
495 : /* return p-adic roots of f, precision prec */
496 : GEN
497 48 : polrootspadic(GEN f, GEN Tp, long prec)
498 : {
499 48 : pari_sp av = avma;
500 : GEN lead, y, T, p;
501 : long PREC, i, k, v;
502 : int reverse;
503 :
504 48 : if (!ff_parse_Tp(Tp, &T,&p,0)) pari_err_TYPE("polrootspadic",Tp);
505 48 : if (typ(f)!=t_POL) pari_err_TYPE("polrootspadic",f);
506 48 : if (gequal0(f)) pari_err_ROOTS0("polrootspadic");
507 48 : if (prec <= 0)
508 6 : pari_err_DOMAIN("polrootspadic", "precision", "<=",gen_0,stoi(prec));
509 42 : f = T? QpXQX_to_ZXY(f, p): QpX_to_ZX(f, p);
510 42 : v = RgX_valrem(f, &f);
511 42 : f = pnormalize(f, p, T, prec, 1, &lead, &PREC, &reverse);
512 42 : y = T? ZXY_ZpQ_roots(f,T,p,PREC): ZX_Zp_roots(f,p,PREC);
513 42 : k = lg(y);
514 42 : if (lead != gen_1) y = RgC_Rg_div(y, lead);
515 42 : if (reverse)
516 6 : for (i=1; i<k; i++) gel(y,i) = ginv(gel(y,i));
517 42 : if (v) y = shallowconcat(zeropadic_shallow(p, prec), y);
518 42 : return gerepilecopy(av, y);
519 : }
520 :
521 : /*******************************************************************/
522 : /* */
523 : /* FACTORIZATION in Zp[X], using ROUND4 */
524 : /* */
525 : /*******************************************************************/
526 :
527 : int
528 2632 : cmp_padic(GEN x, GEN y)
529 : {
530 : long vx, vy;
531 2632 : if (x == gen_0) return -1;
532 2632 : if (y == gen_0) return 1;
533 2632 : vx = valp(x);
534 2632 : vy = valp(y);
535 2632 : if (vx < vy) return 1;
536 2602 : if (vx > vy) return -1;
537 2380 : 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 42 : famat_flatten(GEN fa)
544 : {
545 42 : GEN y, P = gel(fa,1), E = gel(fa,2);
546 42 : long i, l = lg(E);
547 42 : y = cgetg(l, t_VEC);
548 126 : for (i = 1; i < l; i++)
549 : {
550 84 : GEN p = gel(P,i);
551 84 : long e = itou(gel(E,i));
552 84 : gel(y,i) = const_col(e, p);
553 : }
554 42 : y = shallowconcat1(y); return mkmat2(y, const_col(lg(y)-1, gen_1));
555 : }
556 :
557 : GEN
558 498 : factorpadic(GEN f, GEN p, long r)
559 : {
560 498 : pari_sp av = avma;
561 : GEN y, ppow;
562 : long v, n;
563 498 : int reverse = 0, exact;
564 :
565 498 : if (typ(f)!=t_POL) pari_err_TYPE("factorpadic",f);
566 498 : if (typ(p)!=t_INT) pari_err_TYPE("factorpadic",p);
567 498 : if (r <= 0) pari_err_DOMAIN("factorpadic", "precision", "<=",gen_0,stoi(r));
568 492 : if (!signe(f)) return prime_fact(f);
569 492 : if (!degpol(f)) return trivial_fact();
570 :
571 492 : exact = RgX_is_QX(f); /* before RgX_valrem which may lose type information */
572 492 : v = RgX_valrem_inexact(f, &f);
573 492 : ppow = powiu(p,r);
574 492 : n = degpol(f);
575 492 : if (!n)
576 24 : y = trivial_fact();
577 : else
578 : {
579 : GEN P, lead;
580 : long i, l, pr;
581 :
582 468 : f = QpX_to_ZX(f, p);
583 468 : f = pnormalize(f, p, NULL, r, n-1, &lead, &pr, &reverse);
584 468 : y = ZpX_monic_factor(f, p, pr);
585 468 : P = gel(y,1); l = lg(P);
586 468 : if (!isint1(lead))
587 252 : for (i=1; i<l; i++) gel(P,i) = Q_primpart(RgX_unscale(gel(P,i), lead));
588 1578 : for (i=1; i<l; i++)
589 : {
590 1110 : GEN t = gel(P,i);
591 1110 : if (reverse) t = RgX_recip_shallow(t);
592 1110 : gel(P,i) = ZX_to_ZpX_normalized(t,p,ppow,r);
593 : }
594 : }
595 492 : if (v)
596 : { /* v > 0 */
597 54 : GEN X = ZX_to_ZpX(pol_x(varn(f)), p, ppow, r);
598 54 : y = famat_mulpow_shallow(y, X, utoipos(v));
599 : }
600 492 : if (!exact) y = famat_flatten(y);
601 492 : return gerepilecopy(av, sort_factor_pol(y, cmp_padic));
602 : }
|
