Line data Source code
1 : /* Copyright (C) 2000 The PARI group.
2 :
3 : This file is part of the PARI/GP package.
4 :
5 : PARI/GP is free software; you can redistribute it and/or modify it under the
6 : terms of the GNU General Public License as published by the Free Software
7 : Foundation; 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 : #include "pari.h"
16 : #include "paripriv.h"
17 : /*********************************************************************/
18 : /** PERFECT POWERS **/
19 : /*********************************************************************/
20 : #define DEBUGLEVEL DEBUGLEVEL_arith
21 :
22 : /*********************************************************************/
23 : /** INTEGRAL LOGARITHM **/
24 : /*********************************************************************/
25 : /* y > 1 and B > 0 integers. Return e such that y^e <= B < y^(e+1), i.e
26 : * e = floor(log_y B). Set *ptq = y^e if non-NULL */
27 : long
28 833142 : ulogintall(ulong B, ulong y, ulong *ptq)
29 : {
30 : ulong r, r2;
31 : long e;
32 :
33 833142 : if (y == 2)
34 : {
35 23067 : long eB = expu(B); /* 2^eB <= B < 2^(eB + 1) */
36 23067 : if (ptq) *ptq = 1UL << eB;
37 23067 : return eB;
38 : }
39 810075 : r = y, r2 = 1UL;
40 2933611 : for (e=1;; e++)
41 : { /* here, r = y^e, r2 = y^(e-1) */
42 2933611 : if (r >= B)
43 : {
44 808263 : if (r != B) { e--; r = r2; }
45 808263 : if (ptq) *ptq = r;
46 808263 : return e;
47 : }
48 2125348 : r2 = r;
49 2125348 : r = umuluu_or_0(y, r);
50 2125380 : if (!r)
51 : {
52 1844 : if (ptq) *ptq = r2;
53 1844 : return e;
54 : }
55 : }
56 : }
57 :
58 : /* y > 1 and B > 0 integers. Return e such that y^e <= B < y^(e+1), i.e
59 : * e = floor(log_y B). Set *ptq = y^e if non-NULL */
60 : long
61 934044 : logintall(GEN B, GEN y, GEN *ptq)
62 : {
63 : pari_sp av;
64 934044 : long ey, e, emax, i, eB = expi(B); /* 2^eB <= B < 2^(eB + 1) */
65 : GEN q, pow2;
66 :
67 934045 : if (lgefint(B) == 3)
68 : {
69 : ulong q;
70 831993 : if (lgefint(y) > 3)
71 : {
72 0 : if (ptq) *ptq = gen_1;
73 0 : return 0;
74 : }
75 831993 : if (!ptq) return ulogintall(B[2], y[2], NULL);
76 127363 : e = ulogintall(B[2], y[2], &q);
77 127367 : *ptq = utoi(q); return e;
78 : }
79 102052 : if (equaliu(y,2))
80 : {
81 872 : if (ptq) *ptq = int2n(eB);
82 872 : return eB;
83 : }
84 101190 : av = avma;
85 101190 : ey = expi(y);
86 : /* eB/(ey+1) - 1 < e <= eB/ey */
87 101192 : emax = eB/ey;
88 101192 : if (emax <= 13) /* e small, be naive */
89 : {
90 10894 : GEN r = y, r2 = gen_1;
91 10894 : for (e=1;; e++)
92 105335 : { /* here, r = y^e, r2 = y^(e-1) */
93 116229 : long fl = cmpii(r, B);
94 116229 : if (fl >= 0)
95 : {
96 10894 : if (fl) { e--; cgiv(r); r = r2; }
97 10894 : if (ptq) *ptq = gerepileuptoint(av, r); else set_avma(av);
98 10894 : return e;
99 : }
100 105335 : r2 = r; r = mulii(r,y);
101 : }
102 : }
103 : /* e >= 13 ey / (ey+1) >= 6.5 */
104 :
105 : /* binary splitting: compute bits of e one by one */
106 : /* compute pow2[i] = y^(2^i) [i < crude upper bound for log_2 log_y(B)] */
107 90298 : pow2 = new_chunk((long)log2(eB)+2);
108 90299 : gel(pow2,0) = y;
109 90299 : for (i=0, q=y;; )
110 422823 : {
111 513122 : GEN r = gel(pow2,i); /* r = y^2^i */
112 513122 : long fl = cmpii(r,B);
113 513090 : if (!fl)
114 : {
115 0 : e = 1L<<i;
116 0 : if (ptq) *ptq = gerepileuptoint(av, r); else set_avma(av);
117 0 : return e;
118 : }
119 513090 : if (fl > 0) { i--; break; }
120 481461 : q = r;
121 481461 : if (1L<<(i+1) > emax) break;
122 422792 : gel(pow2,++i) = sqri(q);
123 : }
124 :
125 90298 : for (e = 1L<<i;;)
126 391175 : { /* y^e = q < B < r = q * y^(2^i) */
127 481473 : pari_sp av2 = avma;
128 : long fl;
129 : GEN r;
130 481473 : if (--i < 0) break;
131 391187 : r = mulii(q, gel(pow2,i));
132 391189 : fl = cmpii(r, B);
133 391221 : if (fl > 0) set_avma(av2);
134 : else
135 : {
136 176014 : e += (1L<<i);
137 176014 : q = r;
138 176014 : if (!fl) break; /* B = r */
139 : }
140 : }
141 90293 : if (ptq) *ptq = gerepileuptoint(av, q); else set_avma(av);
142 90294 : return e;
143 : }
144 :
145 : long
146 26322 : logint0(GEN B, GEN y, GEN *ptq)
147 : {
148 26322 : const char *f = "logint";
149 26322 : if (typ(y) != t_INT) pari_err_TYPE(f,y);
150 26322 : if (cmpis(y, 2) < 0) pari_err_DOMAIN(f, "b" ,"<=", gen_1, y);
151 26322 : if (typ(B) != t_INT)
152 : {
153 77 : pari_sp av = avma;
154 : long a;
155 77 : if (typ(B) == t_REAL)
156 : {
157 : long e, p;
158 35 : if (cmprs(B, 1) < 1) pari_err_DOMAIN(f, "x", "<", gen_1, B);
159 28 : e = expo(B); if (e < 0) return 0;
160 28 : if (equaliu(y, 2)) return e;
161 14 : if (expu(e) < 50)
162 : {
163 14 : a = floor(dbllog2(B) / dbllog2(y));
164 14 : if (ptq) *ptq = powiu(y, a);
165 14 : return a;
166 : }
167 : /* play safe */
168 0 : p = lg(B);
169 0 : if (nbits2lg(e+1) > p)
170 : { /* try to avoid precision loss in truncation */
171 0 : if (p > DEFAULTPREC) { p = DEFAULTPREC; B = rtor(B, p); }
172 0 : a = itos(floorr(divrr(logr_abs(B), logr_abs(itor(y, p)))));
173 0 : set_avma(av); if (ptq) *ptq = powiu(y, a);
174 0 : return a;
175 : }
176 0 : a = logintall(truncr(B), y, ptq);
177 : }
178 : else
179 : {
180 42 : GEN b = gfloor(B);
181 42 : if (typ(b) != t_INT) pari_err_TYPE(f,B);
182 42 : if (signe(b) <= 0) pari_err_DOMAIN(f, "x", "<", gen_1, B);
183 28 : a = logintall(b, y, ptq);
184 : }
185 28 : if (!ptq) return gc_long(av, a);
186 14 : *ptq = gerepileuptoint(av, *ptq); return a;
187 : }
188 26245 : if (signe(B) <= 0) pari_err_DOMAIN(f, "x" ,"<=", gen_0, B);
189 26245 : return logintall(B,y,ptq);
190 : }
191 :
192 : /*********************************************************************/
193 : /** INTEGRAL SQUARE ROOT **/
194 : /*********************************************************************/
195 : GEN
196 88457 : sqrtint(GEN a)
197 : {
198 88457 : if (typ(a) != t_INT)
199 : {
200 56 : pari_sp av = avma;
201 56 : if (typ(a) == t_REAL)
202 : {
203 : long e;
204 28 : switch(signe(a))
205 : {
206 0 : case 0: return gen_0;
207 7 : case -1: pari_err_DOMAIN("sqrtint", "argument", "<", gen_0,a);
208 : }
209 21 : e = expo(a); if (e < 0) return gen_0;
210 21 : if (nbits2lg(e+1) > lg(a))
211 14 : a = floorr(sqrtr(a)); /* try to avoid precision loss in truncation */
212 : else
213 7 : a = sqrti(truncr(a));
214 : }
215 : else
216 : {
217 28 : GEN b = gfloor(a);
218 28 : if (typ(b) != t_INT) pari_err_TYPE("sqrtint",a);
219 21 : if (signe(b) < 0) pari_err_DOMAIN("sqrtint", "argument", "<", gen_0,a);
220 14 : a = sqrti(b);
221 : }
222 28 : return gerepileuptoleaf(av, a);
223 : }
224 88401 : switch (signe(a))
225 : {
226 88380 : case 1: return sqrti(a);
227 7 : case 0: return gen_0;
228 14 : default: pari_err_DOMAIN("sqrtint", "argument", "<", gen_0,a);
229 : }
230 : return NULL; /* LCOV_EXCL_LINE */
231 : }
232 : GEN
233 119 : sqrtint0(GEN a, GEN *r)
234 : {
235 119 : if (!r) return sqrtint(a);
236 49 : if (typ(a) != t_INT)
237 : {
238 28 : GEN b = sqrtint(a);
239 28 : pari_sp av = avma;
240 28 : *r = gerepileupto(av, gsub(a, sqri(b))); return b;
241 : }
242 21 : switch (signe(a))
243 : {
244 14 : case 1: return sqrtremi(a, r);
245 7 : case 0: *r = gen_0; return gen_0;
246 0 : default: pari_err_DOMAIN("sqrtint", "argument", "<", gen_0,a);
247 : }
248 : return NULL; /* LCOV_EXCL_LINE */
249 : }
250 :
251 : /*********************************************************************/
252 : /** PERFECT SQUARE **/
253 : /*********************************************************************/
254 : static int
255 14399628 : squaremod(ulong A)
256 : {
257 14399628 : const int squaremod64[]={
258 : 1,1,0,0,1,0,0,0,0,1, 0,0,0,0,0,0,1,1,0,0, 0,0,0,0,0,1,0,0,0,0,
259 : 0,0,0,1,0,0,1,0,0,0, 0,1,0,0,0,0,0,0,0,1, 0,0,0,0,0,0,0,1,0,0, 0,0,0,0};
260 14399628 : const int squaremod63[]={
261 : 1,1,0,0,1,0,0,1,0,1, 0,0,0,0,0,0,1,0,1,0, 0,0,1,0,0,1,0,0,1,0,
262 : 0,0,0,0,0,0,1,1,0,0, 0,0,0,1,0,0,1,0,0,1, 0,0,0,0,0,0,0,0,1,0, 0,0,0};
263 14399628 : const int squaremod65[]={
264 : 1,1,0,0,1,0,0,0,0,1, 1,0,0,0,1,0,1,0,0,0, 0,0,0,0,0,1,1,0,0,1,
265 : 1,0,0,0,0,1,1,0,0,1, 1,0,0,0,0,0,0,0,0,1, 0,1,0,0,0,1,1,0,0,0, 0,1,0,0,1};
266 14399628 : const int squaremod11[]={1,1,0,1,1,1,0,0,0,1, 0};
267 14399628 : return (squaremod64[A & 0x3fUL]
268 6777761 : && squaremod63[A % 63UL]
269 4711560 : && squaremod65[A % 65UL]
270 21177389 : && squaremod11[A % 11UL]);
271 : }
272 :
273 : /* emulate Z_issquareall on single-word integers */
274 : long
275 10762270 : uissquareall(ulong A, ulong *sqrtA)
276 : {
277 10762270 : if (!A) { *sqrtA = 0; return 1; }
278 10762270 : if (squaremod(A))
279 : {
280 1937789 : ulong a = usqrt(A);
281 1937788 : if (a * a == A) { *sqrtA = a; return 1; }
282 : }
283 8857534 : return 0;
284 : }
285 : long
286 274804 : uissquare(ulong A)
287 : {
288 274804 : if (!A) return 1;
289 274804 : if (squaremod(A)) { ulong a = usqrt(A); if (a * a == A) return 1; }
290 265008 : return 0;
291 : }
292 :
293 : long
294 6681360 : Z_issquareall(GEN x, GEN *pt)
295 : {
296 : pari_sp av;
297 : GEN y, r;
298 :
299 6681360 : switch(signe(x))
300 : {
301 181009 : case -1: return 0;
302 1813 : case 0: if (pt) *pt=gen_0; return 1;
303 : }
304 6498538 : if (lgefint(x) == 3)
305 : {
306 3135991 : ulong u = uel(x,2), a;
307 3135991 : if (!pt) return uissquare(u);
308 2861187 : if (!uissquareall(u, &a)) return 0;
309 1476560 : *pt = utoipos(a); return 1;
310 : }
311 3362547 : if (!squaremod(umodiu(x, 64*63*65*11))) return 0;
312 2211712 : av = avma; y = sqrtremi(x, &r);
313 2211712 : if (r != gen_0) return gc_long(av,0);
314 41058 : if (pt) { *pt = y; set_avma((pari_sp)y); } else set_avma(av);
315 41058 : return 1;
316 : }
317 :
318 : /* a t_INT, p prime */
319 : long
320 111180 : Zp_issquare(GEN a, GEN p)
321 : {
322 : long v;
323 : GEN ap;
324 :
325 111180 : if (!signe(a) || equali1(a)) return 1;
326 111180 : v = Z_pvalrem(a, p, &ap);
327 111180 : if (v&1) return 0;
328 63425 : return absequaliu(p, 2)? Mod8(ap) == 1
329 63425 : : kronecker(ap,p) == 1;
330 : }
331 :
332 : static long
333 3780 : polissquareall(GEN x, GEN *pt)
334 : {
335 : pari_sp av;
336 : long v;
337 : GEN y, a, b, p;
338 :
339 3780 : if (!signe(x))
340 : {
341 7 : if (pt) *pt = gcopy(x);
342 7 : return 1;
343 : }
344 3773 : if (odd(degpol(x))) return 0; /* odd degree */
345 3773 : av = avma;
346 3773 : v = RgX_valrem(x, &x);
347 3773 : if (v & 1) return gc_long(av,0);
348 3766 : a = gel(x,2); /* test constant coeff */
349 3766 : if (!pt)
350 63 : { if (!issquare(a)) return gc_long(av,0); }
351 : else
352 3703 : { if (!issquareall(a,&b)) return gc_long(av,0); }
353 3766 : if (!degpol(x)) { /* constant polynomial */
354 63 : if (!pt) return gc_long(av,1);
355 28 : y = scalarpol(b, varn(x)); goto END;
356 : }
357 3703 : p = characteristic(x);
358 3703 : if (signe(p) && !mod2(p))
359 : {
360 : long i, lx;
361 35 : if (!absequaliu(p,2)) pari_err_IMPL("issquare for even characteristic != 2");
362 28 : x = gmul(x, mkintmod(gen_1, gen_2));
363 28 : lx = lg(x);
364 28 : if ((lx-3) & 1) return gc_long(av,0);
365 49 : for (i = 3; i < lx; i+=2)
366 28 : if (!gequal0(gel(x,i))) return gc_long(av,0);
367 21 : if (pt) {
368 14 : y = cgetg((lx+3) / 2, t_POL);
369 49 : for (i = 2; i < lx; i+=2)
370 35 : if (!issquareall(gel(x,i), &gel(y, (i+2)>>1))) return gc_long(av,0);
371 14 : y[1] = evalsigne(1) | evalvarn(varn(x));
372 14 : goto END;
373 : } else {
374 21 : for (i = 2; i < lx; i+=2)
375 14 : if (!issquare(gel(x,i))) return gc_long(av,0);
376 7 : return gc_long(av,1);
377 : }
378 : }
379 : else
380 : {
381 3668 : long m = 1;
382 3668 : x = RgX_Rg_div(x,a);
383 : /* a(x^m) = B^2 => B = b(x^m) provided a(0) != 0 */
384 3668 : if (!signe(p)) x = RgX_deflate_max(x,&m);
385 3668 : y = ser2rfrac_i(gsqrt(RgX_to_ser(x,lg(x)-1),0));
386 3675 : if (!RgX_equal(RgX_sqr(y), x)) return gc_long(av,0);
387 3661 : if (!pt) return gc_long(av,1);
388 3654 : if (!gequal1(a)) y = gmul(b, y);
389 3654 : if (m != 1) y = RgX_inflate(y,m);
390 : }
391 3696 : END:
392 3696 : if (v) y = RgX_shift_shallow(y, v>>1);
393 3696 : *pt = gerepilecopy(av, y); return 1;
394 : }
395 :
396 : /* b unit mod p */
397 : static int
398 651 : Up_ispower(GEN b, GEN K, GEN p, long d, GEN *pt)
399 : {
400 651 : if (d == 1)
401 : { /* mod p: faster */
402 567 : if (!Fp_ispower(b, K, p)) return 0;
403 560 : if (pt) *pt = Fp_sqrtn(b, K, p, NULL);
404 : }
405 : else
406 : { /* mod p^{2 +} */
407 84 : if (!ispower(cvtop(b, p, d), K, pt)) return 0;
408 63 : if (pt) *pt = gtrunc(*pt);
409 : }
410 623 : return 1;
411 : }
412 :
413 : /* We're studying whether a mod (q*p^e) is a K-th power, (q,p) = 1.
414 : * Decide mod p^e, then reduce a mod q unless q = NULL. */
415 : static int
416 798 : handle_pe(GEN *pa, GEN q, GEN L, GEN K, GEN p, long e)
417 : {
418 : GEN t, A;
419 798 : long v = Z_pvalrem(*pa, p, &A), d = e - v;
420 798 : if (d <= 0) t = gen_0;
421 : else
422 : {
423 : ulong r;
424 735 : v = uabsdivui_rem(v, K, &r);
425 735 : if (r || !Up_ispower(A, K, p, d, L? &t: NULL)) return 0;
426 623 : if (L && v) t = mulii(t, powiu(p, v));
427 : }
428 686 : if (q) *pa = modii(*pa, q);
429 686 : if (L) vectrunc_append(L, mkintmod(t, powiu(p, e)));
430 686 : return 1;
431 : }
432 : long
433 756 : Zn_ispower(GEN a, GEN q, GEN K, GEN *pt)
434 : {
435 : GEN L, N;
436 : pari_sp av;
437 : long e, i, l;
438 : ulong pp;
439 : forprime_t S;
440 :
441 756 : if (!signe(a))
442 : {
443 91 : if (pt) {
444 91 : GEN t = cgetg(3, t_INTMOD);
445 91 : gel(t,1) = icopy(q); gel(t,2) = gen_0; *pt = t;
446 : }
447 91 : return 1;
448 : }
449 : /* a != 0 */
450 665 : av = avma;
451 :
452 665 : if (typ(q) != t_INT) /* integer factorization */
453 : {
454 0 : GEN P = gel(q,1), E = gel(q,2);
455 0 : l = lg(P);
456 0 : L = pt? vectrunc_init(l): NULL;
457 0 : for (i = 1; i < l; i++)
458 : {
459 0 : GEN p = gel(P,i);
460 0 : long e = itos(gel(E,i));
461 0 : if (!handle_pe(&a, NULL, L, K, p, e)) return gc_long(av,0);
462 : }
463 0 : goto END;
464 : }
465 665 : if (!mod2(K)
466 546 : && kronecker(a, shifti(q,-vali(q))) == -1) return gc_long(av,0);
467 658 : L = pt? vectrunc_init(expi(q)+1): NULL;
468 658 : u_forprime_init(&S, 2, tridiv_bound(q));
469 4128228 : while ((pp = u_forprime_next(&S)))
470 : {
471 : int stop;
472 4128067 : e = Z_lvalrem_stop(&q, pp, &stop);
473 4128067 : if (!e) continue;
474 539 : if (!handle_pe(&a, q, L, K, utoipos(pp), e)) return gc_long(av,0);
475 518 : if (stop)
476 : {
477 476 : if (!is_pm1(q) && !handle_pe(&a, q, L, K, q, 1)) return gc_long(av,0);
478 476 : goto END;
479 : }
480 : }
481 161 : l = lg(primetab);
482 161 : for (i = 1; i < l; i++)
483 : {
484 0 : GEN p = gel(primetab,i);
485 0 : e = Z_pvalrem(q, p, &q);
486 0 : if (!e) continue;
487 0 : if (!handle_pe(&a, q, L, K, p, e)) return gc_long(av,0);
488 0 : if (is_pm1(q)) goto END;
489 : }
490 161 : N = gcdii(a,q);
491 161 : if (!is_pm1(N))
492 : {
493 112 : if (ifac_isprime(N))
494 : {
495 70 : e = Z_pvalrem(q, N, &q);
496 70 : if (!handle_pe(&a, q, L, K, N, e)) return gc_long(av,0);
497 : }
498 : else
499 : {
500 42 : GEN part = ifac_start(N, 0);
501 : for(;;)
502 42 : {
503 : long e;
504 : GEN p;
505 84 : if (!ifac_next(&part, &p, &e)) break;
506 42 : e = Z_pvalrem(q, p, &q);
507 42 : if (!handle_pe(&a, q, L, K, p, e)) return gc_long(av,0);
508 : }
509 : }
510 : }
511 91 : if (!is_pm1(q))
512 : {
513 91 : if (ifac_isprime(q))
514 : {
515 28 : if (!handle_pe(&a, q, L, K, q, 1)) return gc_long(av,0);
516 : }
517 : else
518 : { /* icopy needed: ifac_next would destroy q */
519 63 : GEN part = ifac_start(icopy(q), 0);
520 : for(;;)
521 91 : {
522 : long e;
523 : GEN p;
524 154 : if (!ifac_next(&part, &p, &e)) break;
525 112 : if (!handle_pe(&a, q, L, K, p, e)) return gc_long(av,0);
526 : }
527 : }
528 : }
529 0 : END:
530 546 : if (pt) *pt = gerepileupto(av, chinese1_coprime_Z(L));
531 546 : return 1;
532 : }
533 :
534 : static long
535 56 : polmodispower(GEN x, GEN K, GEN *pt)
536 : {
537 56 : pari_sp av = avma;
538 56 : GEN p = NULL, T = NULL;
539 56 : if (Rg_is_FpXQ(x, &T,&p) && p)
540 : {
541 42 : x = liftall_shallow(x);
542 42 : if (T) T = liftall_shallow(T);
543 42 : if (!Fq_ispower(x, K, T, p)) return gc_long(av,0);
544 28 : if (!pt) return gc_long(av,1);
545 21 : x = Fq_sqrtn(x, K, T,p, NULL);
546 21 : if (typ(x) == t_INT)
547 7 : x = Fp_to_mod(x,p);
548 : else
549 14 : x = mkpolmod(FpX_to_mod(x,p), FpX_to_mod(T,p));
550 21 : *pt = gerepilecopy(av, x); return 1;
551 : }
552 14 : pari_err_IMPL("ispower for general t_POLMOD");
553 0 : return 0;
554 : }
555 : static long
556 56 : rfracispower(GEN x, GEN K, GEN *pt)
557 : {
558 56 : pari_sp av = avma;
559 56 : GEN n = gel(x,1), d = gel(x,2);
560 56 : long v = -RgX_valrem(d, &d), vx = varn(d);
561 56 : if (typ(n) == t_POL && varn(n) == vx) v += RgX_valrem(n, &n);
562 56 : if (!dvdsi(v, K)) return gc_long(av, 0);
563 49 : if (lg(d) >= 3)
564 : {
565 49 : GEN a = gel(d,2); /* constant term */
566 49 : if (!gequal1(a)) { d = RgX_Rg_div(d, a); n = gdiv(n, a); }
567 : }
568 49 : if (!ispower(d, K, pt? &d: NULL) || !ispower(n, K, pt? &n: NULL))
569 0 : return gc_long(av, 0);
570 49 : if (!pt) return gc_long(av, 1);
571 28 : x = gdiv(n, d);
572 28 : if (v) x = gmul(x, monomial(gen_1, v / itos(K), vx));
573 28 : *pt = gerepileupto(av, x); return 1;
574 : }
575 : long
576 252268 : issquareall(GEN x, GEN *pt)
577 : {
578 252268 : long tx = typ(x);
579 : GEN F;
580 : pari_sp av;
581 :
582 252268 : if (!pt) return issquare(x);
583 38086 : switch(tx)
584 : {
585 17499 : case t_INT: return Z_issquareall(x, pt);
586 1939 : case t_FRAC: av = avma;
587 1939 : F = cgetg(3, t_FRAC);
588 1939 : if ( !Z_issquareall(gel(x,1), &gel(F,1))
589 1939 : || !Z_issquareall(gel(x,2), &gel(F,2))) return gc_long(av,0);
590 1876 : *pt = F; return 1;
591 :
592 21 : case t_POLMOD:
593 21 : return polmodispower(x, gen_2, pt);
594 3710 : case t_POL: return polissquareall(x,pt);
595 21 : case t_RFRAC: return rfracispower(x, gen_2, pt);
596 :
597 14791 : case t_REAL: case t_COMPLEX: case t_PADIC: case t_SER:
598 14791 : if (!issquare(x)) return 0;
599 14791 : *pt = gsqrt(x, DEFAULTPREC); return 1;
600 :
601 63 : case t_INTMOD:
602 63 : return Zn_ispower(gel(x,2), gel(x,1), gen_2, pt);
603 :
604 42 : case t_FFELT: return FF_issquareall(x, pt);
605 :
606 : }
607 0 : pari_err_TYPE("issquareall",x);
608 : return 0; /* LCOV_EXCL_LINE */
609 : }
610 :
611 : long
612 229267 : issquare(GEN x)
613 : {
614 : GEN a, p;
615 : long v;
616 :
617 229267 : switch(typ(x))
618 : {
619 214077 : case t_INT:
620 214077 : return Z_issquare(x);
621 :
622 14721 : case t_REAL:
623 14721 : return (signe(x)>=0);
624 :
625 84 : case t_INTMOD:
626 84 : return Zn_ispower(gel(x,2), gel(x,1), gen_2, NULL);
627 :
628 35 : case t_FRAC:
629 35 : return Z_issquare(gel(x,1)) && Z_issquare(gel(x,2));
630 :
631 7 : case t_FFELT: return FF_issquareall(x, NULL);
632 :
633 56 : case t_COMPLEX:
634 56 : return 1;
635 :
636 133 : case t_PADIC:
637 133 : a = gel(x,4); if (!signe(a)) return 1;
638 133 : if (valp(x)&1) return 0;
639 119 : p = gel(x,2);
640 119 : if (!absequaliu(p, 2)) return (kronecker(a,p) != -1);
641 :
642 42 : v = precp(x); /* here p=2, a is odd */
643 42 : if ((v>=3 && mod8(a) != 1 ) ||
644 21 : (v==2 && mod4(a) != 1)) return 0;
645 21 : return 1;
646 :
647 21 : case t_POLMOD:
648 21 : return polmodispower(x, gen_2, NULL);
649 :
650 70 : case t_POL:
651 70 : return polissquareall(x,NULL);
652 :
653 49 : case t_SER:
654 49 : if (!signe(x)) return 1;
655 42 : if (valser(x)&1) return 0;
656 35 : return issquare(gel(x,2));
657 :
658 14 : case t_RFRAC:
659 14 : return rfracispower(x, gen_2, NULL);
660 : }
661 0 : pari_err_TYPE("issquare",x);
662 : return 0; /* LCOV_EXCL_LINE */
663 : }
664 0 : GEN gissquare(GEN x) { return issquare(x)? gen_1: gen_0; }
665 0 : GEN gissquareall(GEN x, GEN *pt) { return issquareall(x,pt)? gen_1: gen_0; }
666 :
667 : long
668 1386 : ispolygonal(GEN x, GEN S, GEN *N)
669 : {
670 1386 : pari_sp av = avma;
671 : GEN D, d, n;
672 1386 : if (typ(x) != t_INT) pari_err_TYPE("ispolygonal", x);
673 1386 : if (typ(S) != t_INT) pari_err_TYPE("ispolygonal", S);
674 1386 : if (abscmpiu(S,3) < 0) pari_err_DOMAIN("ispolygonal","s","<", utoipos(3),S);
675 1386 : if (signe(x) < 0) return 0;
676 1386 : if (signe(x) == 0) { if (N) *N = gen_0; return 1; }
677 1260 : if (is_pm1(x)) { if (N) *N = gen_1; return 1; }
678 : /* n = (sqrt( (8s - 16) x + (s-4)^2 ) + s - 4) / 2(s - 2) */
679 1134 : if (abscmpiu(S, 1<<16) < 0) /* common case ! */
680 : {
681 441 : ulong s = S[2], r;
682 441 : if (s == 4) return Z_issquareall(x, N);
683 378 : if (s == 3)
684 0 : D = addiu(shifti(x, 3), 1);
685 : else
686 378 : D = addiu(mului(8*s - 16, x), (s-4)*(s-4));
687 378 : if (!Z_issquareall(D, &d)) return gc_long(av,0);
688 378 : if (s == 3)
689 0 : d = subiu(d, 1);
690 : else
691 378 : d = addiu(d, s - 4);
692 378 : n = absdiviu_rem(d, 2*s - 4, &r);
693 378 : if (r) return gc_long(av,0);
694 : }
695 : else
696 : {
697 693 : GEN r, S_2 = subiu(S,2), S_4 = subiu(S,4);
698 693 : D = addii(mulii(shifti(S_2,3), x), sqri(S_4));
699 693 : if (!Z_issquareall(D, &d)) return gc_long(av,0);
700 693 : d = addii(d, S_4);
701 693 : n = dvmdii(shifti(d,-1), S_2, &r);
702 693 : if (r != gen_0) return gc_long(av,0);
703 : }
704 1071 : if (N) *N = gerepileuptoint(av, n); else set_avma(av);
705 1071 : return 1;
706 : }
707 :
708 : /*********************************************************************/
709 : /** PERFECT POWER **/
710 : /*********************************************************************/
711 : static long
712 952 : polispower(GEN x, GEN K, GEN *pt)
713 : {
714 : pari_sp av;
715 952 : long v, d, k = itos(K);
716 : GEN y, a, b;
717 952 : GEN T = NULL, p = NULL;
718 :
719 952 : if (!signe(x))
720 : {
721 7 : if (pt) *pt = gcopy(x);
722 7 : return 1;
723 : }
724 945 : d = degpol(x);
725 945 : if (d % k) return 0; /* degree not multiple of k */
726 938 : av = avma;
727 938 : if (RgX_is_FpXQX(x, &T, &p) && p)
728 : { /* over Fq */
729 336 : if (T && typ(T) == t_FFELT)
730 : {
731 126 : if (!FFX_ispower(x, k, T, pt)) return gc_long(av,0);
732 105 : return 1;
733 : }
734 210 : x = RgX_to_FqX(x,T,p);
735 210 : if (!FqX_ispower(x, k, T,p, pt)) return gc_long(av,0);
736 175 : if (pt) *pt = gerepileupto(av, FqX_to_mod(*pt, T, p));
737 175 : return 1;
738 : }
739 602 : v = RgX_valrem(x, &x);
740 602 : if (v % k) return 0;
741 595 : v /= k;
742 595 : a = gel(x,2); b = NULL;
743 595 : if (!ispower(a, K, &b)) return gc_long(av,0);
744 581 : if (d)
745 : {
746 553 : GEN p = characteristic(x);
747 553 : a = leading_coeff(x);
748 553 : if (!ispower(a, K, &b)) return gc_long(av,0);
749 553 : x = RgX_normalize(x);
750 553 : if (signe(p) && cmpii(p,K) <= 0)
751 0 : pari_err_IMPL("ispower(general t_POL) in small characteristic");
752 553 : y = gtrunc(gsqrtn(RgX_to_ser(x,lg(x)), K, NULL, 0));
753 553 : if (!RgX_equal(powgi(y, K), x)) return gc_long(av,0);
754 : }
755 : else
756 28 : y = pol_1(varn(x));
757 581 : if (pt)
758 : {
759 546 : if (!gequal1(a))
760 : {
761 35 : if (!b) b = gsqrtn(a, K, NULL, DEFAULTPREC);
762 35 : y = gmul(b,y);
763 : }
764 546 : if (v) y = RgX_shift_shallow(y, v);
765 546 : *pt = gerepilecopy(av, y);
766 : }
767 35 : else set_avma(av);
768 581 : return 1;
769 : }
770 :
771 : long
772 1559237 : Z_ispowerall(GEN x, ulong k, GEN *pt)
773 : {
774 1559237 : long s = signe(x);
775 : ulong mask;
776 1559237 : if (!s) { if (pt) *pt = gen_0; return 1; }
777 1559237 : if (s > 0) {
778 1558908 : if (k == 2) return Z_issquareall(x, pt);
779 1425582 : if (k == 3) { mask = 1; return !!is_357_power(x, pt, &mask); }
780 227878 : if (k == 5) { mask = 2; return !!is_357_power(x, pt, &mask); }
781 60305 : if (k == 7) { mask = 4; return !!is_357_power(x, pt, &mask); }
782 54187 : return is_kth_power(x, k, pt);
783 : }
784 329 : if (!odd(k)) return 0;
785 140 : if (Z_ispowerall(absi_shallow(x), k, pt))
786 : {
787 126 : if (pt) *pt = negi(*pt);
788 126 : return 1;
789 : };
790 14 : return 0;
791 : }
792 :
793 : /* is x a K-th power mod p ? Assume p prime. */
794 : int
795 567 : Fp_ispower(GEN x, GEN K, GEN p)
796 : {
797 567 : pari_sp av = avma;
798 : GEN p_1;
799 567 : x = modii(x, p);
800 567 : if (!signe(x) || equali1(x)) return gc_bool(av,1);
801 : /* implies p > 2 */
802 126 : p_1 = subiu(p,1);
803 126 : K = gcdii(K, p_1);
804 126 : if (absequaliu(K, 2)) return gc_bool(av, kronecker(x,p) > 0);
805 49 : x = Fp_pow(x, diviiexact(p_1,K), p);
806 49 : return gc_bool(av, equali1(x));
807 : }
808 :
809 : /* x unit defined modulo 2^e, e > 0, p prime */
810 : static int
811 2541 : U2_issquare(GEN x, long e)
812 : {
813 2541 : long r = signe(x)>=0?mod8(x):8-mod8(x);
814 2541 : if (e==1) return 1;
815 2541 : if (e==2) return (r&3L) == 1;
816 2009 : return r == 1;
817 : }
818 : /* x unit defined modulo p^e, e > 0, p prime */
819 : static int
820 5229 : Up_issquare(GEN x, GEN p, long e)
821 5229 : { return (absequaliu(p,2))? U2_issquare(x, e): kronecker(x,p)==1; }
822 :
823 : long
824 2793 : Zn_issquare(GEN d, GEN fn)
825 : {
826 : long j, np;
827 2793 : if (typ(d) != t_INT) pari_err_TYPE("Zn_issquare",d);
828 2793 : if (typ(fn) == t_INT) return Zn_ispower(d, fn, gen_2, NULL);
829 : /* integer factorization */
830 2793 : np = nbrows(fn);
831 6006 : for (j = 1; j <= np; ++j)
832 : {
833 5586 : GEN r, p = gcoeff(fn, j, 1);
834 5586 : long e = itos(gcoeff(fn, j, 2));
835 5586 : long v = Z_pvalrem(d,p,&r);
836 5586 : if (v < e && (odd(v) || !Up_issquare(r, p, e-v))) return 0;
837 : }
838 420 : return 1;
839 : }
840 :
841 : static long
842 1113 : Qp_ispower(GEN x, GEN K, GEN *pt)
843 : {
844 1113 : pari_sp av = avma;
845 1113 : GEN z = Qp_sqrtn(x, K, NULL);
846 1113 : if (!z) return gc_long(av,0);
847 819 : if (pt) *pt = z;
848 819 : return 1;
849 : }
850 :
851 : long
852 8555983 : ispower(GEN x, GEN K, GEN *pt)
853 : {
854 : GEN z;
855 :
856 8555983 : if (!K) return gisanypower(x, pt);
857 1555801 : if (typ(K) != t_INT) pari_err_TYPE("ispower",K);
858 1555801 : if (signe(K) <= 0) pari_err_DOMAIN("ispower","exponent","<=",gen_0,K);
859 1555801 : if (equali1(K)) { if (pt) *pt = gcopy(x); return 1; }
860 1555738 : switch(typ(x)) {
861 1482796 : case t_INT:
862 1482796 : if (lgefint(K) != 3) return 0;
863 1482740 : return Z_ispowerall(x, itou(K), pt);
864 70184 : case t_FRAC:
865 : {
866 70184 : GEN a = gel(x,1), b = gel(x,2);
867 : ulong k;
868 70184 : if (lgefint(K) != 3) return 0;
869 70177 : k = itou(K);
870 70177 : if (pt) {
871 70114 : z = cgetg(3, t_FRAC);
872 70114 : if (Z_ispowerall(a, k, &a) && Z_ispowerall(b, k, &b)) {
873 1484 : *pt = z; gel(z,1) = a; gel(z,2) = b; return 1;
874 : }
875 68630 : set_avma((pari_sp)(z + 3)); return 0;
876 : }
877 63 : return Z_ispower(a, k) && Z_ispower(b, k);
878 : }
879 609 : case t_INTMOD:
880 609 : return Zn_ispower(gel(x,2), gel(x,1), K, pt);
881 28 : case t_FFELT:
882 28 : return FF_ispower(x, K, pt);
883 :
884 1113 : case t_PADIC:
885 1113 : return Qp_ispower(x, K, pt);
886 14 : case t_POLMOD:
887 14 : return polmodispower(x, K, pt);
888 952 : case t_POL:
889 952 : return polispower(x, K, pt);
890 21 : case t_RFRAC:
891 21 : return rfracispower(x, K, pt);
892 7 : case t_REAL:
893 7 : if (signe(x) < 0 && !mpodd(K)) return 0;
894 : case t_COMPLEX:
895 14 : if (pt) *pt = gsqrtn(x, K, NULL, DEFAULTPREC);
896 14 : return 1;
897 :
898 7 : case t_SER:
899 7 : if (signe(x) && (!dvdsi(valser(x), K) || !ispower(gel(x,2), K, NULL)))
900 0 : return 0;
901 7 : if (pt) *pt = gsqrtn(x, K, NULL, DEFAULTPREC);
902 7 : return 1;
903 : }
904 0 : pari_err_TYPE("ispower",x);
905 : return 0; /* LCOV_EXCL_LINE */
906 : }
907 :
908 : long
909 7000182 : gisanypower(GEN x, GEN *pty)
910 : {
911 7000182 : long tx = typ(x);
912 : ulong k, h;
913 7000182 : if (tx == t_INT) return Z_isanypower(x, pty);
914 14 : if (tx == t_FRAC)
915 : {
916 14 : pari_sp av = avma;
917 14 : GEN fa, P, E, a = gel(x,1), b = gel(x,2);
918 : long i, j, p, e;
919 14 : int sw = (abscmpii(a, b) > 0);
920 :
921 14 : if (sw) swap(a, b);
922 14 : k = Z_isanypower(a, pty? &a: NULL);
923 14 : if (!k)
924 : { /* a = -1,1 or not a pure power */
925 7 : if (!is_pm1(a)) return gc_long(av,0);
926 7 : if (signe(a) < 0) b = negi(b);
927 7 : k = Z_isanypower(b, pty? &b: NULL);
928 7 : if (!k || !pty) return gc_long(av,k);
929 7 : *pty = gerepileupto(av, ginv(b));
930 7 : return k;
931 : }
932 7 : fa = factoru(k);
933 7 : P = gel(fa,1);
934 7 : E = gel(fa,2); h = k;
935 14 : for (i = lg(P) - 1; i > 0; i--)
936 : {
937 7 : p = P[i];
938 7 : e = E[i];
939 21 : for (j = 0; j < e; j++)
940 14 : if (!is_kth_power(b, p, &b)) break;
941 7 : if (j < e) k /= upowuu(p, e - j);
942 : }
943 7 : if (k == 1) return gc_long(av,0);
944 7 : if (!pty) return gc_long(av,k);
945 0 : if (k != h) a = powiu(a, h/k);
946 0 : *pty = gerepilecopy(av, mkfrac(a, b));
947 0 : return k;
948 : }
949 0 : pari_err_TYPE("gisanypower", x);
950 : return 0; /* LCOV_EXCL_LINE */
951 : }
952 :
953 : /* v_p(x) = e != 0 for some p; return ispower(x,,&x), updating x.
954 : * No need to optimize for 2,3,5,7 powers (done before) */
955 : static long
956 507443 : split_exponent(ulong e, GEN *x)
957 : {
958 : GEN fa, P, E;
959 507443 : long i, j, l, k = 1;
960 507443 : if (e == 1) return 1;
961 14 : fa = factoru(e);
962 14 : P = gel(fa,1);
963 14 : E = gel(fa,2); l = lg(P);
964 28 : for (i = 1; i < l; i++)
965 : {
966 14 : ulong p = P[i];
967 28 : for (j = 0; j < E[i]; j++)
968 : {
969 : GEN y;
970 14 : if (!is_kth_power(*x, p, &y)) break;
971 14 : k *= p; *x = y;
972 : }
973 : }
974 14 : return k;
975 : }
976 :
977 : /* any prime divisor of x is > 102 */
978 : static long
979 869912 : Z_isanypower_101(GEN *px)
980 : {
981 869912 : const double LOG2_103 = 6.6865; /* lower bound for log_2(103) */
982 869912 : const double LOG103 = 4.6347; /* lower bound for log(103) */
983 : forprime_t T;
984 869912 : ulong mask = 7, e2;
985 : long k, ex;
986 869912 : GEN y, x = *px;
987 :
988 869912 : k = 1;
989 873099 : while (Z_issquareall(x, &y)) { k <<= 1; x = y; }
990 870340 : while ( (ex = is_357_power(x, &y, &mask)) ) { k *= ex; x = y; }
991 869912 : e2 = (ulong)((expi(x) + 1) / LOG2_103); /* >= log_103 (x) */
992 869912 : if (u_forprime_init(&T, 11, e2))
993 : {
994 17367 : GEN logx = NULL;
995 17367 : const ulong Q = 30011; /* prime */
996 : ulong p, xmodQ;
997 17367 : double dlogx = 0;
998 : /* cut off at x^(1/p) ~ 2^30 bits which seems to be about optimum;
999 : * for large p the modular checks are no longer competitively fast */
1000 17409 : while ( (ex = is_pth_power(x, &y, &T, 30)) )
1001 : {
1002 42 : k *= ex; x = y;
1003 42 : e2 = (ulong)((expi(x) + 1) / LOG2_103);
1004 42 : u_forprime_restrict(&T, e2);
1005 : }
1006 17367 : if (DEBUGLEVEL>4)
1007 0 : err_printf("Z_isanypower: now k=%ld, x=%ld-bit\n", k, expi(x)+1);
1008 17367 : xmodQ = umodiu(x, Q);
1009 : /* test Q | x, just in case */
1010 17367 : if (!xmodQ) { *px = x; return k * split_exponent(Z_lval(x,Q), px); }
1011 : /* x^(1/p) < 2^31 */
1012 17353 : p = T.p;
1013 17353 : if (p <= e2)
1014 : {
1015 17339 : logx = logr_abs( itor(x, DEFAULTPREC) );
1016 17339 : dlogx = rtodbl(logx);
1017 17339 : e2 = (ulong)(dlogx / LOG103); /* >= log_103(x) */
1018 : }
1019 143129 : while (p && p <= e2)
1020 : { /* is x a p-th power ? By computing y = round(x^(1/p)).
1021 : * Check whether y^p = x, first mod Q, then exactly. */
1022 125776 : pari_sp av = avma;
1023 : long e;
1024 125776 : GEN logy = divru(logx, p), y = grndtoi(mpexp(logy), &e);
1025 125776 : ulong ymodQ = umodiu(y,Q);
1026 125776 : if (e >= -10 || Fl_powu(ymodQ, p % (Q-1), Q) != xmodQ
1027 125776 : || !equalii(powiu(y, p), x)) set_avma(av);
1028 : else
1029 : {
1030 21 : k *= p; x = y; xmodQ = ymodQ; logx = logy; dlogx /= p;
1031 21 : e2 = (ulong)(dlogx / LOG103); /* >= log_103(x) */
1032 21 : u_forprime_restrict(&T, e2);
1033 21 : continue; /* if success, retry same p */
1034 : }
1035 125755 : p = u_forprime_next(&T);
1036 : }
1037 : }
1038 869898 : *px = x; return k;
1039 : }
1040 :
1041 : static ulong tinyprimes[] = {
1042 : 2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71,
1043 : 73, 79, 83, 89, 97, 101, 103, 107, 109, 113, 127, 131, 137, 139, 149, 151,
1044 : 157, 163, 167, 173, 179, 181, 191, 193, 197, 199
1045 : };
1046 :
1047 : /* disregard the sign of x, caller will take care of x < 0 */
1048 : static long
1049 7195587 : Z_isanypower_aux(GEN x, GEN *pty)
1050 : {
1051 : long ex, v, i, l, k;
1052 : GEN y, P, E;
1053 7195587 : ulong mask, e = 0;
1054 :
1055 7195587 : if (abscmpii(x, gen_2) < 0) return 0; /* -1,0,1 */
1056 :
1057 7194849 : x = absi_shallow(x);
1058 7194849 : k = l = 1;
1059 7194849 : P = cgetg(26 + 1, t_VECSMALL);
1060 7194849 : E = cgetg(26 + 1, t_VECSMALL);
1061 : /* trial division */
1062 61980961 : for(i = 0; i < 26; i++)
1063 : {
1064 60639502 : ulong p = tinyprimes[i];
1065 : int stop;
1066 60639502 : v = Z_lvalrem_stop(&x, p, &stop);
1067 60639502 : if (v)
1068 : {
1069 8043961 : P[l] = p;
1070 8043961 : E[l] = v; l++;
1071 8424015 : e = ugcd(e, v); if (e == 1) goto END;
1072 : }
1073 55166166 : if (stop) {
1074 380054 : if (is_pm1(x)) k = e;
1075 380054 : goto END;
1076 : }
1077 : }
1078 :
1079 1341459 : if (e)
1080 : { /* Bingo. Result divides e */
1081 : long v3, v5, v7;
1082 507429 : ulong e2 = e;
1083 507429 : v = u_lvalrem(e2, 2, &e2);
1084 507429 : if (v)
1085 : {
1086 377780 : for (i = 0; i < v; i++)
1087 : {
1088 375659 : if (!Z_issquareall(x, &y)) break;
1089 2333 : k <<= 1; x = y;
1090 : }
1091 : }
1092 507429 : mask = 0;
1093 507429 : v3 = u_lvalrem(e2, 3, &e2); if (v3) mask = 1;
1094 507429 : v5 = u_lvalrem(e2, 5, &e2); if (v5) mask |= 2;
1095 507429 : v7 = u_lvalrem(e2, 7, &e2); if (v7) mask |= 4;
1096 507506 : while ( (ex = is_357_power(x, &y, &mask)) ) {
1097 77 : x = y;
1098 77 : switch(ex)
1099 : {
1100 28 : case 3: k *= 3; if (--v3 == 0) mask &= ~1; break;
1101 28 : case 5: k *= 5; if (--v5 == 0) mask &= ~2; break;
1102 21 : case 7: k *= 7; if (--v7 == 0) mask &= ~4; break;
1103 : }
1104 : }
1105 507429 : k *= split_exponent(e2, &x);
1106 : }
1107 : else
1108 834030 : k = Z_isanypower_101(&x);
1109 7194849 : END:
1110 7194849 : if (pty && k != 1)
1111 : {
1112 69933 : if (e)
1113 : { /* add missing small factors */
1114 66865 : y = powuu(P[1], E[1] / k);
1115 76860 : for (i = 2; i < l; i++) y = mulii(y, powuu(P[i], E[i] / k));
1116 66865 : x = equali1(x)? y: mulii(x,y);
1117 : }
1118 69933 : *pty = x;
1119 : }
1120 7194849 : return k == 1? 0: k;
1121 : }
1122 :
1123 : long
1124 7195587 : Z_isanypower(GEN x, GEN *pty)
1125 : {
1126 7195587 : pari_sp av = avma;
1127 7195587 : long k = Z_isanypower_aux(x, pty);
1128 7195587 : if (!k) return gc_long(av,0);
1129 69996 : if (signe(x) < 0)
1130 : {
1131 42 : long v = vals(k);
1132 42 : if (v)
1133 : {
1134 28 : k >>= v;
1135 28 : if (k == 1) return gc_long(av,0);
1136 21 : if (!pty) return gc_long(av,k);
1137 14 : *pty = gerepileuptoint(av, powiu(*pty, 1<<v));
1138 14 : togglesign(*pty); return k;
1139 : }
1140 14 : if (pty) togglesign_safe(pty);
1141 : }
1142 69968 : if (!pty) return gc_long(av, k);
1143 69912 : *pty = gerepilecopy(av, *pty); return k;
1144 : }
1145 :
1146 : /* Faster than expi(n) == vali(n) or hamming(n) == 1 even for single-word
1147 : * values. If all you have is a word, you can just use n & !(n & (n-1)). */
1148 : long
1149 528555 : Z_ispow2(GEN n)
1150 : {
1151 : GEN xp;
1152 : long i, l;
1153 : ulong u;
1154 528555 : if (signe(n) != 1) return 0;
1155 472548 : xp = int_LSW(n); u = *xp; l = lgefint(n);
1156 533429 : for (i = 3; i < l; ++i)
1157 : {
1158 252507 : if (u) return 0;
1159 60881 : xp = int_nextW(xp); u = *xp;
1160 : }
1161 280922 : return !(u & (u-1));
1162 : }
1163 :
1164 : static long
1165 842139 : isprimepower_i(GEN n, GEN *pt, long flag)
1166 : {
1167 842139 : pari_sp av = avma;
1168 : long i, v;
1169 :
1170 842139 : if (typ(n) != t_INT) pari_err_TYPE("isprimepower", n);
1171 842139 : if (signe(n) <= 0) return 0;
1172 :
1173 842139 : if (lgefint(n) == 3)
1174 : {
1175 : ulong p;
1176 541232 : v = uisprimepower(n[2], &p);
1177 541232 : if (v)
1178 : {
1179 55006 : if (pt) *pt = utoipos(p);
1180 55006 : return v;
1181 : }
1182 486226 : return 0;
1183 : }
1184 1663684 : for (i = 0; i < 26; i++)
1185 : {
1186 1627802 : ulong p = tinyprimes[i];
1187 1627802 : v = Z_lvalrem(n, p, &n);
1188 1627802 : if (v)
1189 : {
1190 265025 : set_avma(av);
1191 265025 : if (!is_pm1(n)) return 0;
1192 629 : if (pt) *pt = utoipos(p);
1193 629 : return v;
1194 : }
1195 : }
1196 : /* p | n => p >= 103 */
1197 35882 : v = Z_isanypower_101(&n); /* expensive */
1198 35882 : if (!(flag? isprime(n): BPSW_psp(n))) return gc_long(av,0);
1199 5570 : if (pt) *pt = gerepilecopy(av, n); else set_avma(av);
1200 5570 : return v;
1201 : }
1202 : long
1203 840147 : isprimepower(GEN n, GEN *pt) { return isprimepower_i(n,pt,1); }
1204 : long
1205 1992 : ispseudoprimepower(GEN n, GEN *pt) { return isprimepower_i(n,pt,0); }
1206 :
1207 : long
1208 642865 : uisprimepower(ulong n, ulong *pp)
1209 : { /* We must have CUTOFF^11 >= ULONG_MAX and CUTOFF^3 < ULONG_MAX.
1210 : * Tests suggest that 200-300 is the best range for 64-bit platforms. */
1211 642865 : const ulong CUTOFF = 200UL;
1212 642865 : const long TINYCUTOFF = 46; /* tinyprimes[45] = 199 */
1213 642865 : const ulong CUTOFF3 = CUTOFF*CUTOFF*CUTOFF;
1214 : #ifdef LONG_IS_64BIT
1215 : /* primes preceeding the appropriate root of ULONG_MAX. */
1216 568182 : const ulong ROOT9 = 137;
1217 568182 : const ulong ROOT8 = 251;
1218 568182 : const ulong ROOT7 = 563;
1219 568182 : const ulong ROOT5 = 7129;
1220 568182 : const ulong ROOT4 = 65521;
1221 : #else
1222 74683 : const ulong ROOT9 = 11;
1223 74683 : const ulong ROOT8 = 13;
1224 74683 : const ulong ROOT7 = 23;
1225 74683 : const ulong ROOT5 = 83;
1226 74683 : const ulong ROOT4 = 251;
1227 : #endif
1228 : ulong mask;
1229 : long v, i;
1230 : int e;
1231 642865 : if (n < 2) return 0;
1232 642823 : if (!odd(n)) {
1233 341774 : if (n & (n-1)) return 0;
1234 65474 : *pp = 2; return vals(n);
1235 : }
1236 301048 : if (n < 8) { *pp = n; return 1; } /* 3,5,7 */
1237 3655110 : for (i = 1/*skip p=2*/; i < TINYCUTOFF; i++)
1238 : {
1239 3596035 : ulong p = tinyprimes[i];
1240 3596035 : if (n % p == 0)
1241 : {
1242 212447 : v = u_lvalrem(n, p, &n);
1243 212447 : if (n == 1) { *pp = p; return v; }
1244 210012 : return 0;
1245 : }
1246 : }
1247 : /* p | n => p >= CUTOFF */
1248 :
1249 59075 : if (n < CUTOFF3)
1250 : {
1251 46354 : if (n < CUTOFF*CUTOFF || uisprime_101(n)) { *pp = n; return 1; }
1252 0 : if (uissquareall(n, &n)) { *pp = n; return 2; }
1253 0 : return 0;
1254 : }
1255 :
1256 : /* Check for squares, fourth powers, and eighth powers as appropriate. */
1257 12721 : v = 1;
1258 12721 : if (uissquareall(n, &n)) {
1259 0 : v <<= 1;
1260 0 : if (CUTOFF <= ROOT4 && uissquareall(n, &n)) {
1261 0 : v <<= 1;
1262 0 : if (CUTOFF <= ROOT8 && uissquareall(n, &n)) v <<= 1;
1263 : }
1264 : }
1265 :
1266 12721 : if (CUTOFF > ROOT5) mask = 1;
1267 : else
1268 : {
1269 12720 : const ulong CUTOFF5 = CUTOFF3*CUTOFF*CUTOFF;
1270 12720 : if (n < CUTOFF5) mask = 1; else mask = 3;
1271 12720 : if (CUTOFF <= ROOT7)
1272 : {
1273 12720 : const ulong CUTOFF7 = CUTOFF5*CUTOFF*CUTOFF;
1274 12720 : if (n >= CUTOFF7) mask = 7;
1275 : }
1276 : }
1277 :
1278 12721 : if (CUTOFF <= ROOT9 && (e = uis_357_power(n, &n, &mask))) { v *= e; mask=1; }
1279 12721 : if ((e = uis_357_power(n, &n, &mask))) v *= e;
1280 :
1281 12721 : if (uisprime_101(n)) { *pp = n; return v; }
1282 6984 : return 0;
1283 : }
1284 :
|
