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 12887 : polmod_mod(GEN x, GEN y)
31 : {
32 12887 : GEN z, a, T = gel(x,1);
33 12887 : 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 1054558 : rem_scal_pol(GEN x, GEN y)
59 : {
60 1054558 : if (degpol(y))
61 : {
62 1052990 : if (!signe(y)) pari_err_INV("grem",y);
63 1052990 : 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 4199262 : grem(GEN x, GEN y)
122 : {
123 4199262 : const char *f = "euclidean division";
124 4199262 : long tx = typ(x), ty = typ(y), vx = gvar(x), vy = gvar(y);
125 :
126 4199263 : if (ty == t_POL)
127 : {
128 4199200 : if (varncmp(vx,vy) >= 0)
129 : {
130 : pari_sp av;
131 : GEN z;
132 4199202 : if (!signe(y)) pari_err_INV("grem",y);
133 4199203 : if (vx != vy) return rem_scal_pol(x,y);
134 3144645 : switch(tx)
135 : {
136 12887 : case t_POLMOD: return polmod_mod(x,y);
137 3119627 : 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 = gel(x,2);
188 336 : if (p && !equalii(p, q)) pari_err_MODULUS("Zp_to_Z", p,q);
189 315 : }
190 : /* shallow */
191 : static GEN
192 4326 : Zp_to_Z(GEN x, GEN p) {
193 4326 : switch(typ(x))
194 : {
195 4088 : 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 4305 : return x;
202 : }
203 : /* shallow */
204 : static GEN
205 749 : ZpX_to_ZX(GEN f, GEN p) {
206 749 : long i, l = lg(f);
207 749 : GEN F = cgetg_copy(f, &l); F[1] = f[1];
208 4914 : for (i=2; i<l; i++) gel(F,i) = Zp_to_Z(gel(f,i), p);
209 735 : return F;
210 : }
211 :
212 : static GEN
213 686 : get_padic_content(GEN f, GEN p)
214 : {
215 686 : GEN c = content(f);
216 686 : if (gequal0(c)) /* O(p^n) can occur */
217 : {
218 0 : if (typ(c) != t_PADIC) pari_err_TYPE("QpX_to_ZX",f);
219 0 : check_padic_p(c, p);
220 0 : c = powis(p, valp(c));
221 : }
222 686 : return c;
223 : }
224 : /* make f suitable for [root|factor]padic. Shallow */
225 : static GEN
226 623 : QpX_to_ZX(GEN f, GEN p)
227 : {
228 623 : GEN c = get_padic_content(f, p);
229 623 : f = RgX_Rg_div(f, c);
230 623 : return ZpX_to_ZX(f, p);
231 : }
232 :
233 : /* x in Z return x + O(pr), pr = p^r. Shallow */
234 : static GEN
235 4361 : Z_to_Zp(GEN x, GEN p, GEN pr, long r)
236 : {
237 : GEN y;
238 4361 : long v, sx = signe(x);
239 :
240 4361 : if (!sx) return zeropadic_shallow(p,r);
241 3829 : v = Z_pvalrem(x,p,&x);
242 3829 : if (v) {
243 840 : if (r <= v) return zeropadic_shallow(p,r);
244 735 : r -= v;
245 735 : pr = powiu(p,r);
246 : }
247 3724 : y = cgetg(5,t_PADIC);
248 3724 : y[1] = evalprecp(r)|evalvalp(v);
249 3724 : gel(y,2) = p;
250 3724 : gel(y,3) = pr;
251 3724 : gel(y,4) = modii(x,pr); return y;
252 : }
253 : /* shallow */
254 : static GEN
255 56 : ZV_to_ZpV(GEN z, GEN p, long r)
256 : {
257 56 : long i, l = lg(z);
258 56 : GEN Z = cgetg(l, typ(z)), q = powiu(p, r);
259 161 : for (i=1; i<l; i++) gel(Z,i) = Z_to_Zp(gel(z,i),p,q,r);
260 56 : return Z;
261 : }
262 : /* shallow */
263 : static GEN
264 1253 : ZX_to_ZpX(GEN z, GEN p, GEN q, long r)
265 : {
266 1253 : long i, l = lg(z);
267 1253 : GEN Z = cgetg(l, t_POL); Z[1] = z[1];
268 5327 : for (i=2; i<l; i++) gel(Z,i) = Z_to_Zp(gel(z,i),p,q,r);
269 1253 : return Z;
270 : }
271 : /* return (x + O(p^r)) normalized (multiply by a unit such that leading coeff
272 : * is a power of p), x in Z[X] (or Z_p[X]). Shallow */
273 : static GEN
274 1169 : ZX_to_ZpX_normalized(GEN x, GEN p, GEN pr, long r)
275 : {
276 1169 : long i, lx = lg(x);
277 1169 : GEN z, lead = leading_coeff(x);
278 :
279 1169 : if (gequal1(lead)) return ZX_to_ZpX(x, p, pr, r);
280 56 : (void)Z_pvalrem(lead, p, &lead); lead = Fp_inv(lead, pr);
281 56 : z = cgetg(lx,t_POL);
282 238 : for (i=2; i < lx; i++) gel(z,i) = Z_to_Zp(mulii(lead,gel(x,i)),p,pr,r);
283 56 : z[1] = x[1]; return z;
284 : }
285 : static GEN
286 49 : ZXV_to_ZpXQV(GEN z, GEN T, GEN p, long r)
287 : {
288 49 : long i, l = lg(z);
289 49 : GEN Z = cgetg(l, typ(z)), q = powiu(p, r);
290 49 : T = ZX_copy(T);
291 126 : for (i=1; i<lg(z); i++) gel(Z,i) = mkpolmod(ZX_to_ZpX(gel(z,i),p,q,r),T);
292 49 : return Z;
293 : }
294 : /* shallow */
295 : static GEN
296 63 : QpXQX_to_ZXY(GEN f, GEN p)
297 : {
298 63 : GEN c = get_padic_content(f, p);
299 63 : long i, l = lg(f);
300 63 : f = RgX_Rg_div(f,c);
301 287 : for (i=2; i<l; i++)
302 : {
303 231 : GEN t = gel(f,i);
304 231 : switch(typ(t))
305 : {
306 91 : case t_POLMOD:
307 91 : t = gel(t,2);
308 91 : t = (typ(t) == t_POL)? ZpX_to_ZX(t, p): Zp_to_Z(t, p);
309 91 : break;
310 0 : case t_POL: t = ZpX_to_ZX(t, p); break;
311 140 : default: t = Zp_to_Z(t, p); break;
312 : }
313 224 : gel(f,i) = t;
314 : }
315 56 : return f;
316 : }
317 :
318 : /*******************************************************************/
319 : /* */
320 : /* p-ADIC ROOTS */
321 : /* */
322 : /*******************************************************************/
323 :
324 : /* f primitive ZX, squarefree, leading term prime to p; 0 <= a < p such that
325 : * f(a) = 0 mod p. Return p-adic roots of f equal to a mod p, in
326 : * precision >= prec */
327 : GEN
328 2863 : ZX_Zp_root(GEN f, GEN a, GEN p, long prec)
329 : {
330 : GEN z, R;
331 : long i, j, k;
332 :
333 2863 : if (signe(FpX_eval(FpX_deriv(f, p), a, p)))
334 : { /* simple zero mod p, go all the way to p^prec */
335 2632 : if (prec > 1) a = ZpX_liftroot(f, a, p, prec);
336 2632 : return mkcol(a);
337 : }
338 :
339 231 : f = ZX_unscale_div(ZX_translate(f,a), p); /* f(pX + a) / p */
340 231 : (void)ZX_pvalrem(f,p,&f);
341 231 : z = cgetg(degpol(f)+1,t_COL);
342 :
343 231 : R = FpX_roots(f, p);
344 574 : for (j=i=1; i<lg(R); i++)
345 : {
346 343 : GEN u = ZX_Zp_root(f, gel(R,i), p, prec-1);
347 756 : for (k=1; k<lg(u); k++) gel(z,j++) = addii(a, mulii(p, gel(u,k)));
348 : }
349 231 : setlg(z,j); return z;
350 : }
351 :
352 : /* a t_PADIC, return vector of p-adic roots of f equal to a (mod p) */
353 : GEN
354 56 : Zp_appr(GEN f, GEN a)
355 : {
356 56 : pari_sp av = avma;
357 56 : GEN z, p = gel(a,2);
358 56 : long v = valp(a), prec = v;
359 :
360 56 : if (signe(gel(a,4))) prec += precp(a);
361 56 : f = QpX_to_ZX(f, p);
362 42 : if (degpol(f) <= 0) pari_err_CONSTPOL("Zp_appr");
363 42 : if (v < 0) pari_err_DOMAIN("padicappr", "v(a)", "<", gen_0, stoi(v));
364 35 : f = ZX_radical(f);
365 35 : a = padic_to_Fp(a, p);
366 35 : if (signe(FpX_eval(f,a,p))) { set_avma(av); return cgetg(1,t_COL); }
367 28 : z = ZX_Zp_root(f, a, p, prec);
368 28 : return gerepilecopy(av, ZV_to_ZpV(z, p, prec));
369 : }
370 : static long
371 126 : pval(GEN x, GEN p) { return typ(x) == t_INT? Z_pval(x,p): ZX_pval(x,p); }
372 : /* f a ZXX, f(0) != 0 */
373 : static GEN
374 539 : pnormalize(GEN f, GEN p, GEN T, long prec, long n,
375 : GEN *plead, long *pprec, int *prev)
376 : {
377 539 : *plead = leading_coeff(f);
378 539 : *pprec = prec;
379 539 : *prev = 0;
380 539 : if (!isint1(*plead))
381 : {
382 63 : long v = pval(*plead,p), v1 = pval(constant_coeff(f),p);
383 63 : if (v1 < v)
384 : {
385 42 : *prev = 1;
386 42 : f = RgX_recip_i(f); /* f(0) != 0 so degree is the same */
387 : /* beware loss of precision from lc(factor), whose valuation is <= v */
388 42 : *pprec += v; v = v1;
389 : }
390 63 : *pprec += v * n;
391 : }
392 539 : if (!T) return ZX_Q_normalize(f, plead);
393 14 : *plead = gen_1;
394 14 : return FpXQX_normalize(f, T, powiu(p,*pprec));
395 : }
396 :
397 : /**************************************************************************/
398 :
399 : static void
400 238 : scalar_getprec(GEN x, long *pprec, GEN *pp)
401 : {
402 238 : if (typ(x)==t_PADIC)
403 : {
404 98 : long e = valp(x); if (signe(gel(x,4))) e += precp(x);
405 98 : if (e < *pprec) *pprec = e;
406 98 : check_padic_p(x, *pp);
407 98 : *pp = gel(x,2);
408 : }
409 238 : }
410 : static void
411 98 : getprec(GEN x, long *pprec, GEN *pp)
412 : {
413 : long i;
414 98 : if (typ(x) != t_POL) scalar_getprec(x, pprec, pp);
415 : else
416 266 : for (i = lg(x)-1; i>1; i--) scalar_getprec(gel(x,i), pprec, pp);
417 98 : }
418 :
419 : /* assume f(a) = 0 (mod T,p) */
420 : static GEN
421 105 : ZXY_ZpQ_root(GEN f, GEN a, GEN T, GEN p, long prec)
422 : {
423 : GEN z, R;
424 : long i, j, k, lR;
425 105 : if (signe(FqX_eval(FqX_deriv(f,T,p), a, T,p)))
426 : { /* simple zero mod (T,p), go all the way to p^prec */
427 77 : if (prec > 1) a = ZpXQX_liftroot(f, a, T, p, prec);
428 77 : return mkcol(a);
429 : }
430 28 : f = RgX_unscale(RgXQX_translate(f, a, T), p);
431 28 : f = RgX_Rg_div(f, powiu(p, gvaluation(f,p)));
432 28 : z = cgetg(degpol(f)+1,t_COL);
433 28 : R = FpXQX_roots(FqX_red(f,T,p), T, p); lR = lg(R);
434 70 : for(j=i=1; i<lR; i++)
435 : {
436 42 : GEN u = ZXY_ZpQ_root(f, gel(R,i), T, p, prec-1);
437 84 : for (k=1; k<lg(u); k++) gel(z,j++) = gadd(a, gmul(p, gel(u,k)));
438 : }
439 28 : setlg(z,j); return z;
440 : }
441 :
442 : /* a belongs to finite extension of Q_p, return all roots of f equal to a
443 : * mod p. Don't assume f(a) = 0 (mod p) */
444 : GEN
445 105 : padicappr(GEN f, GEN a)
446 : {
447 : GEN p, z, T, Tp;
448 : long prec;
449 105 : pari_sp av = avma;
450 :
451 105 : if (typ(f)!=t_POL) pari_err_TYPE("padicappr",f);
452 105 : switch(typ(a)) {
453 56 : case t_PADIC: return Zp_appr(f,a);
454 49 : case t_POLMOD: break;
455 0 : default: pari_err_TYPE("padicappr",a);
456 : }
457 49 : if (gequal0(f)) pari_err_ROOTS0("padicappr");
458 49 : T = gel(a,1);
459 49 : a = gel(a,2);
460 49 : p = NULL; prec = LONG_MAX;
461 49 : getprec(a, &prec, &p);
462 49 : getprec(T, &prec, &p); if (!p) pari_err_TYPE("padicappr",T);
463 49 : f = QpXQX_to_ZXY(f, p);
464 42 : if (typ(a) != t_POL) a = scalarpol_shallow(a, varn(T));
465 42 : a = ZpX_to_ZX(a,p);
466 42 : T = QpX_to_ZX(T,p);
467 : /* ensure that f /= (f,f') is separable */
468 42 : (void)nfgcd_all(f, RgX_deriv(f), T, NULL, &f);
469 :
470 42 : Tp = FpX_red(T, p); a = FqX_red(a, Tp, p);
471 42 : if (!gequal0(FqX_eval(FqX_red(f,Tp,p), a, Tp,p)))
472 7 : { set_avma(av); return cgetg(1,t_COL); } /* f(a) != 0 (mod p,T) */
473 35 : z = ZXY_ZpQ_root(f, a, T, p, prec);
474 35 : return gerepilecopy(av, ZXV_to_ZpXQV(z, T, p, prec));
475 : }
476 :
477 : /* vector of p-adic roots of the ZX f, leading term prime to p. Shallow */
478 : static GEN
479 35 : ZX_Zp_roots(GEN f, GEN p, long prec)
480 : {
481 : long l, i;
482 : GEN r;
483 :
484 35 : f = ZX_radical(f);
485 35 : r = FpX_roots(f, p);
486 35 : l = lg(r); if (l == 1) return r;
487 91 : for (i = 1; i < l; i++) gel(r,i) = ZX_Zp_root(f, gel(r,i), p, prec);
488 28 : return ZV_to_ZpV(shallowconcat1(r), p, prec);
489 : }
490 : /* vector of p-adic roots of the ZXX f, leading term prime to p. Shallow */
491 : static GEN
492 14 : ZXY_ZpQ_roots(GEN f, GEN T, GEN p, long prec)
493 : {
494 : long l, i;
495 : GEN r;
496 :
497 14 : (void)nfgcd_all(f, RgX_deriv(f), T, NULL, &f);
498 14 : r = FqX_roots(f, FpX_red(T,p), p);
499 14 : l = lg(r); if (l == 1) return r;
500 42 : for (i = 1; i < l; i++) gel(r,i) = ZXY_ZpQ_root(f, gel(r,i), T, p, prec);
501 14 : return ZXV_to_ZpXQV(shallowconcat1(r), T, p, prec);
502 : }
503 :
504 : /* return p-adic roots of f, precision prec */
505 : GEN
506 56 : polrootspadic(GEN f, GEN Tp, long prec)
507 : {
508 56 : pari_sp av = avma;
509 : GEN lead, y, T, p;
510 : long PREC, i, k, v;
511 : int reverse;
512 :
513 56 : if (!ff_parse_Tp(Tp, &T,&p,0)) pari_err_TYPE("polrootspadic",Tp);
514 56 : if (typ(f)!=t_POL) pari_err_TYPE("polrootspadic",f);
515 56 : if (gequal0(f)) pari_err_ROOTS0("polrootspadic");
516 56 : if (prec <= 0)
517 7 : pari_err_DOMAIN("polrootspadic", "precision", "<=",gen_0,stoi(prec));
518 49 : f = T? QpXQX_to_ZXY(f, p): QpX_to_ZX(f, p);
519 49 : v = RgX_valrem(f, &f);
520 49 : f = pnormalize(f, p, T, prec, 1, &lead, &PREC, &reverse);
521 49 : y = T? ZXY_ZpQ_roots(f,T,p,PREC): ZX_Zp_roots(f,p,PREC);
522 49 : k = lg(y);
523 49 : if (lead != gen_1) y = RgC_Rg_div(y, lead);
524 49 : if (reverse)
525 7 : for (i=1; i<k; i++) gel(y,i) = ginv(gel(y,i));
526 49 : if (v) y = shallowconcat(zeropadic_shallow(p, prec), y);
527 49 : return gerepilecopy(av, y);
528 : }
529 :
530 : /*******************************************************************/
531 : /* */
532 : /* FACTORIZATION in Zp[X], using ROUND4 */
533 : /* */
534 : /*******************************************************************/
535 :
536 : int
537 3007 : cmp_padic(GEN x, GEN y)
538 : {
539 : long vx, vy;
540 3007 : if (x == gen_0) return -1;
541 3007 : if (y == gen_0) return 1;
542 3007 : vx = valp(x);
543 3007 : vy = valp(y);
544 3007 : if (vx < vy) return 1;
545 2972 : if (vx > vy) return -1;
546 2734 : return cmpii(gel(x,4), gel(y,4));
547 : }
548 :
549 : /* replace p^e by p*...*p [ factors are not known to be equal, only close at
550 : * input accuracy ] */
551 : static GEN
552 49 : famat_flatten(GEN fa)
553 : {
554 49 : GEN y, P = gel(fa,1), E = gel(fa,2);
555 49 : long i, l = lg(E);
556 49 : y = cgetg(l, t_VEC);
557 147 : for (i = 1; i < l; i++)
558 : {
559 98 : GEN p = gel(P,i);
560 98 : long e = itou(gel(E,i));
561 98 : gel(y,i) = const_col(e, p);
562 : }
563 49 : y = shallowconcat1(y); return mkmat2(y, const_col(lg(y)-1, gen_1));
564 : }
565 :
566 : GEN
567 525 : factorpadic(GEN f, GEN p, long r)
568 : {
569 525 : pari_sp av = avma;
570 : GEN y, ppow;
571 : long v, n;
572 525 : int reverse = 0, exact;
573 :
574 525 : if (typ(f)!=t_POL) pari_err_TYPE("factorpadic",f);
575 525 : if (typ(p)!=t_INT) pari_err_TYPE("factorpadic",p);
576 525 : if (r <= 0) pari_err_DOMAIN("factorpadic", "precision", "<=",gen_0,stoi(r));
577 518 : if (!signe(f)) return prime_fact(f);
578 518 : if (!degpol(f)) return trivial_fact();
579 :
580 518 : exact = RgX_is_QX(f); /* before RgX_valrem which may lose type information */
581 518 : v = RgX_valrem_inexact(f, &f);
582 518 : ppow = powiu(p,r);
583 518 : n = degpol(f);
584 518 : if (!n)
585 28 : y = trivial_fact();
586 : else
587 : {
588 : GEN P, lead;
589 : long i, l, pr;
590 :
591 490 : f = QpX_to_ZX(f, p);
592 490 : f = pnormalize(f, p, NULL, r, n-1, &lead, &pr, &reverse);
593 490 : y = ZpX_monic_factor(f, p, pr);
594 490 : P = gel(y,1); l = lg(P);
595 490 : if (!isint1(lead))
596 266 : for (i=1; i<l; i++) gel(P,i) = Q_primpart(RgX_unscale(gel(P,i), lead));
597 1659 : for (i=1; i<l; i++)
598 : {
599 1169 : GEN t = gel(P,i);
600 1169 : if (reverse) t = RgX_recip_shallow(t);
601 1169 : gel(P,i) = ZX_to_ZpX_normalized(t,p,ppow,r);
602 : }
603 : }
604 518 : if (v)
605 : { /* v > 0 */
606 63 : GEN X = ZX_to_ZpX(pol_x(varn(f)), p, ppow, r);
607 63 : y = famat_mulpow_shallow(y, X, utoipos(v));
608 : }
609 518 : if (!exact) y = famat_flatten(y);
610 518 : return gerepilecopy(av, sort_factor_pol(y, cmp_padic));
611 : }
|
