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 841277 : ulogintall(ulong B, ulong y, ulong *ptq)
29 : {
30 : ulong r, r2;
31 : long e;
32 :
33 841277 : if (y == 2)
34 : {
35 23507 : long eB = expu(B); /* 2^eB <= B < 2^(eB + 1) */
36 23507 : if (ptq) *ptq = 1UL << eB;
37 23507 : return eB;
38 : }
39 817770 : r = y, r2 = 1UL;
40 2963161 : for (e=1;; e++)
41 : { /* here, r = y^e, r2 = y^(e-1) */
42 2963161 : if (r >= B)
43 : {
44 815960 : if (r != B) { e--; r = r2; }
45 815960 : if (ptq) *ptq = r;
46 815960 : return e;
47 : }
48 2147201 : r2 = r;
49 2147201 : r = umuluu_or_0(y, r);
50 2147201 : if (!r)
51 : {
52 1810 : if (ptq) *ptq = r2;
53 1810 : 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 943016 : logintall(GEN B, GEN y, GEN *ptq)
62 : {
63 : pari_sp av;
64 943016 : long ey, e, emax, i, eB = expi(B); /* 2^eB <= B < 2^(eB + 1) */
65 : GEN q, pow2;
66 :
67 943016 : if (lgefint(B) == 3)
68 : {
69 : ulong q;
70 840129 : if (lgefint(y) > 3)
71 : {
72 0 : if (ptq) *ptq = gen_1;
73 0 : return 0;
74 : }
75 840129 : if (!ptq) return ulogintall(B[2], y[2], NULL);
76 130212 : e = ulogintall(B[2], y[2], &q);
77 130212 : *ptq = utoi(q); return e;
78 : }
79 102887 : if (equaliu(y,2))
80 : {
81 941 : if (ptq) *ptq = int2n(eB);
82 941 : return eB;
83 : }
84 101946 : av = avma;
85 101946 : ey = expi(y);
86 : /* eB/(ey+1) - 1 < e <= eB/ey */
87 101946 : emax = eB/ey;
88 101946 : if (emax <= 13) /* e small, be naive */
89 : {
90 11007 : GEN r = y, r2 = gen_1;
91 11007 : for (e=1;; e++)
92 106319 : { /* here, r = y^e, r2 = y^(e-1) */
93 117326 : long fl = cmpii(r, B);
94 117326 : if (fl >= 0)
95 : {
96 11007 : if (fl) { e--; cgiv(r); r = r2; }
97 11007 : if (ptq) *ptq = gc_INT(av, r); else set_avma(av);
98 11007 : return e;
99 : }
100 106319 : 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 90939 : pow2 = new_chunk((long)log2(eB)+2);
108 90939 : gel(pow2,0) = y;
109 90939 : for (i=0, q=y;; )
110 425770 : {
111 516709 : GEN r = gel(pow2,i); /* r = y^2^i */
112 516709 : long fl = cmpii(r,B);
113 516709 : if (!fl)
114 : {
115 0 : e = 1L<<i;
116 0 : if (ptq) *ptq = gc_INT(av, r); else set_avma(av);
117 0 : return e;
118 : }
119 516709 : if (fl > 0) { i--; break; }
120 484778 : q = r;
121 484778 : if (1L<<(i+1) > emax) break;
122 425770 : gel(pow2,++i) = sqri(q);
123 : }
124 :
125 90939 : for (e = 1L<<i;;)
126 393818 : { /* y^e = q < B < r = q * y^(2^i) */
127 484757 : pari_sp av2 = avma;
128 : long fl;
129 : GEN r;
130 484757 : if (--i < 0) break;
131 393825 : r = mulii(q, gel(pow2,i));
132 393825 : fl = cmpii(r, B);
133 393825 : if (fl > 0) set_avma(av2);
134 : else
135 : {
136 177307 : e += (1L<<i);
137 177307 : q = r;
138 177307 : if (!fl) break; /* B = r */
139 : }
140 : }
141 90939 : if (ptq) *ptq = gc_INT(av, q); else set_avma(av);
142 90939 : return e;
143 : }
144 :
145 : long
146 26882 : logint0(GEN B, GEN y, GEN *ptq)
147 : {
148 26882 : const char *f = "logint";
149 26882 : if (typ(y) != t_INT) pari_err_TYPE(f,y);
150 26882 : if (cmpis(y, 2) < 0) pari_err_DOMAIN(f, "b", "<=", gen_1, y);
151 26882 : 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 = gc_INT(av, *ptq); return a;
187 : }
188 26805 : if (signe(B) <= 0) pari_err_DOMAIN(f, "x", "<=", gen_0, B);
189 26805 : return logintall(B,y,ptq);
190 : }
191 :
192 : /*********************************************************************/
193 : /** INTEGRAL SQUARE ROOT **/
194 : /*********************************************************************/
195 : GEN
196 92399 : sqrtint(GEN a)
197 : {
198 92399 : 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 gc_leaf(av, a);
223 : }
224 92343 : switch (signe(a))
225 : {
226 92322 : 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 = gc_upto(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 15862701 : squaremod(ulong A)
256 : {
257 15862701 : 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 15862701 : 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 15862701 : 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 15862701 : const int squaremod11[]={1,1,0,1,1,1,0,0,0,1, 0};
267 15862701 : return (squaremod64[A & 0x3fUL]
268 8041033 : && squaremod63[A % 63UL]
269 5921768 : && squaremod65[A % 65UL]
270 23903734 : && squaremod11[A % 11UL]);
271 : }
272 :
273 : /* emulate Z_issquareall on single-word integers */
274 : long
275 12054475 : uissquareall(ulong A, ulong *sqrtA)
276 : {
277 12054475 : if (!A) { *sqrtA = 0; return 1; }
278 12054475 : if (squaremod(A))
279 : {
280 3135097 : ulong a = usqrt(A);
281 3135097 : if (a * a == A) { *sqrtA = a; return 1; }
282 : }
283 8951362 : return 0;
284 : }
285 : long
286 275003 : uissquare(ulong A)
287 : {
288 275003 : if (!A) return 1;
289 275003 : if (squaremod(A)) { ulong a = usqrt(A); if (a * a == A) return 1; }
290 265204 : return 0;
291 : }
292 :
293 : long
294 8091753 : Z_issquareall(GEN x, GEN *pt)
295 : {
296 : pari_sp av;
297 : GEN y, r;
298 :
299 8091753 : switch(signe(x))
300 : {
301 236512 : case -1: return 0;
302 1631 : case 0: if (pt) *pt=gen_0; return 1;
303 : }
304 7853610 : if (lgefint(x) == 3)
305 : {
306 4320387 : ulong u = uel(x,2), a;
307 4320387 : if (!pt) return uissquare(u);
308 4045384 : if (!uissquareall(u, &a)) return 0;
309 2657724 : *pt = utoipos(a); return 1;
310 : }
311 3533223 : if (!squaremod(umodiu(x, 64*63*65*11))) return 0;
312 2213983 : av = avma; y = sqrtremi(x, &r);
313 2213983 : if (r != gen_0) return gc_long(av,0);
314 40445 : if (pt) { *pt = y; set_avma((pari_sp)y); } else set_avma(av);
315 40445 : return 1;
316 : }
317 :
318 : /* a t_INT, p prime */
319 : long
320 111229 : Zp_issquare(GEN a, GEN p)
321 : {
322 : long v;
323 : GEN ap;
324 :
325 111229 : if (!signe(a) || equali1(a)) return 1;
326 111229 : v = Z_pvalrem(a, p, &ap);
327 111229 : 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 = gc_GEN(av, y); return 1;
394 : }
395 :
396 : static long
397 56 : polmodispower(GEN x, GEN K, GEN *pt)
398 : {
399 56 : pari_sp av = avma;
400 56 : GEN p = NULL, T = NULL;
401 56 : if (Rg_is_FpXQ(x, &T,&p) && p)
402 : {
403 42 : x = liftall_shallow(x);
404 42 : if (T) T = liftall_shallow(T);
405 42 : if (!Fq_ispower(x, K, T, p)) return gc_long(av,0);
406 28 : if (!pt) return gc_long(av,1);
407 21 : x = Fq_sqrtn(x, K, T,p, NULL);
408 21 : if (typ(x) == t_INT)
409 7 : x = Fp_to_mod(x,p);
410 : else
411 14 : x = mkpolmod(FpX_to_mod(x,p), FpX_to_mod(T,p));
412 21 : *pt = gc_GEN(av, x); return 1;
413 : }
414 14 : pari_err_IMPL("ispower for general t_POLMOD");
415 0 : return 0;
416 : }
417 : static long
418 56 : rfracispower(GEN x, GEN K, GEN *pt)
419 : {
420 56 : pari_sp av = avma;
421 56 : GEN n = gel(x,1), d = gel(x,2);
422 56 : long v = -RgX_valrem(d, &d), vx = varn(d);
423 56 : if (typ(n) == t_POL && varn(n) == vx) v += RgX_valrem(n, &n);
424 56 : if (!dvdsi(v, K)) return gc_long(av, 0);
425 49 : if (lg(d) >= 3)
426 : {
427 49 : GEN a = gel(d,2); /* constant term */
428 49 : if (!gequal1(a)) { d = RgX_Rg_div(d, a); n = gdiv(n, a); }
429 : }
430 49 : if (!ispower(d, K, pt? &d: NULL) || !ispower(n, K, pt? &n: NULL))
431 0 : return gc_long(av, 0);
432 49 : if (!pt) return gc_long(av, 1);
433 28 : x = gdiv(n, d);
434 28 : if (v) x = gmul(x, monomial(gen_1, v / itos(K), vx));
435 28 : *pt = gc_upto(av, x); return 1;
436 : }
437 : long
438 433589 : issquareall(GEN x, GEN *pt)
439 : {
440 433589 : long tx = typ(x);
441 : GEN F;
442 : pari_sp av;
443 :
444 433589 : if (!pt) return issquare(x);
445 219407 : switch(tx)
446 : {
447 18164 : case t_INT: return Z_issquareall(x, pt);
448 182595 : case t_FRAC: av = avma;
449 182595 : F = cgetg(3, t_FRAC);
450 182595 : if ( !Z_issquareall(gel(x,1), &gel(F,1))
451 182595 : || !Z_issquareall(gel(x,2), &gel(F,2))) return gc_long(av,0);
452 2275 : *pt = F; return 1;
453 :
454 21 : case t_POLMOD:
455 21 : return polmodispower(x, gen_2, pt);
456 3710 : case t_POL: return polissquareall(x,pt);
457 21 : case t_RFRAC: return rfracispower(x, gen_2, pt);
458 :
459 14791 : case t_REAL: case t_COMPLEX: case t_PADIC: case t_SER:
460 14791 : if (!issquare(x)) return 0;
461 14791 : *pt = gsqrt(x, DEFAULTPREC); return 1;
462 :
463 63 : case t_INTMOD:
464 63 : return Zn_ispower(gel(x,2), gel(x,1), gen_2, pt);
465 :
466 42 : case t_FFELT: return FF_issquareall(x, pt);
467 :
468 : }
469 0 : pari_err_TYPE("issquareall",x);
470 : return 0; /* LCOV_EXCL_LINE */
471 : }
472 :
473 : long
474 229267 : issquare(GEN x)
475 : {
476 : GEN a, p;
477 : long v;
478 :
479 229267 : switch(typ(x))
480 : {
481 214077 : case t_INT:
482 214077 : return Z_issquare(x);
483 :
484 14721 : case t_REAL:
485 14721 : return (signe(x)>=0);
486 :
487 84 : case t_INTMOD:
488 84 : return Zn_issquare(gel(x,2), gel(x,1));
489 :
490 35 : case t_FRAC:
491 35 : return Z_issquare(gel(x,1)) && Z_issquare(gel(x,2));
492 :
493 7 : case t_FFELT: return FF_issquareall(x, NULL);
494 :
495 56 : case t_COMPLEX:
496 56 : return 1;
497 :
498 133 : case t_PADIC:
499 133 : a = padic_u(x); if (!signe(a)) return 1;
500 133 : if (valp(x)&1) return 0;
501 119 : p = padic_p(x);
502 119 : if (!absequaliu(p, 2)) return (kronecker(a,p) != -1);
503 :
504 42 : v = precp(x); /* here p=2, a is odd */
505 42 : if ((v>=3 && mod8(a) != 1 ) ||
506 21 : (v==2 && mod4(a) != 1)) return 0;
507 21 : return 1;
508 :
509 21 : case t_POLMOD:
510 21 : return polmodispower(x, gen_2, NULL);
511 :
512 70 : case t_POL:
513 70 : return polissquareall(x,NULL);
514 :
515 49 : case t_SER:
516 49 : if (!signe(x)) return 1;
517 42 : if (valser(x)&1) return 0;
518 35 : return issquare(gel(x,2));
519 :
520 14 : case t_RFRAC:
521 14 : return rfracispower(x, gen_2, NULL);
522 : }
523 0 : pari_err_TYPE("issquare",x);
524 : return 0; /* LCOV_EXCL_LINE */
525 : }
526 : GEN
527 0 : gissquare(GEN x) { return issquare(x)? gen_1: gen_0; }
528 : GEN
529 0 : gissquareall(GEN x, GEN *pt) { return issquareall(x,pt)? gen_1: gen_0; }
530 :
531 : long
532 1386 : ispolygonal(GEN x, GEN S, GEN *N)
533 : {
534 1386 : pari_sp av = avma;
535 : GEN D, d, n;
536 1386 : if (typ(x) != t_INT) pari_err_TYPE("ispolygonal", x);
537 1386 : if (typ(S) != t_INT) pari_err_TYPE("ispolygonal", S);
538 1386 : if (abscmpiu(S,3) < 0) pari_err_DOMAIN("ispolygonal","s","<", utoipos(3),S);
539 1386 : if (signe(x) < 0) return 0;
540 1386 : if (signe(x) == 0) { if (N) *N = gen_0; return 1; }
541 1260 : if (is_pm1(x)) { if (N) *N = gen_1; return 1; }
542 : /* n = (sqrt( (8s - 16) x + (s-4)^2 ) + s - 4) / 2(s - 2) */
543 1134 : if (abscmpiu(S, 1<<16) < 0) /* common case ! */
544 : {
545 441 : ulong s = S[2], r;
546 441 : if (s == 4) return Z_issquareall(x, N);
547 378 : if (s == 3)
548 0 : D = addiu(shifti(x, 3), 1);
549 : else
550 378 : D = addiu(mului(8*s - 16, x), (s-4)*(s-4));
551 378 : if (!Z_issquareall(D, &d)) return gc_long(av,0);
552 378 : if (s == 3)
553 0 : d = subiu(d, 1);
554 : else
555 378 : d = addiu(d, s - 4);
556 378 : n = absdiviu_rem(d, 2*s - 4, &r);
557 378 : if (r) return gc_long(av,0);
558 : }
559 : else
560 : {
561 693 : GEN r, S_2 = subiu(S,2), S_4 = subiu(S,4);
562 693 : D = addii(mulii(shifti(S_2,3), x), sqri(S_4));
563 693 : if (!Z_issquareall(D, &d)) return gc_long(av,0);
564 693 : d = addii(d, S_4);
565 693 : n = dvmdii(shifti(d,-1), S_2, &r);
566 693 : if (r != gen_0) return gc_long(av,0);
567 : }
568 1071 : if (N) *N = gc_INT(av, n); else set_avma(av);
569 1071 : return 1;
570 : }
571 :
572 : /*********************************************************************/
573 : /** PERFECT POWER **/
574 : /*********************************************************************/
575 : static long
576 1008 : polispower(GEN x, GEN K, GEN *pt)
577 : {
578 : pari_sp av;
579 1008 : long v, d, k = itos(K);
580 : GEN y, a, b;
581 1008 : GEN T = NULL, p = NULL;
582 :
583 1008 : if (!signe(x))
584 : {
585 7 : if (pt) *pt = gcopy(x);
586 7 : return 1;
587 : }
588 1001 : d = degpol(x);
589 1001 : if (d % k) return 0; /* degree not multiple of k */
590 994 : av = avma;
591 994 : if (RgX_is_FpXQX(x, &T, &p) && p)
592 : { /* over Fq */
593 336 : if (T && typ(T) == t_FFELT)
594 : {
595 126 : if (!FFX_ispower(x, k, T, pt)) return gc_long(av,0);
596 105 : return 1;
597 : }
598 210 : x = RgX_to_FqX(x,T,p);
599 210 : if (!FqX_ispower(x, k, T,p, pt)) return gc_long(av,0);
600 175 : if (pt) *pt = gc_upto(av, FqX_to_mod(*pt, T, p));
601 175 : return 1;
602 : }
603 658 : v = RgX_valrem(x, &x);
604 658 : if (v % k) return 0;
605 651 : v /= k;
606 651 : a = gel(x,2); b = NULL;
607 651 : if (!ispower(a, K, &b)) return gc_long(av,0);
608 637 : if (d)
609 : {
610 609 : GEN p = characteristic(x);
611 609 : a = leading_coeff(x);
612 609 : if (!ispower(a, K, &b)) return gc_long(av,0);
613 609 : x = RgX_normalize(x);
614 609 : if (signe(p) && cmpii(p,K) <= 0)
615 0 : pari_err_IMPL("ispower(general t_POL) in small characteristic");
616 609 : y = gtrunc(gsqrtn(RgX_to_ser(x,lg(x)), K, NULL, 0));
617 609 : if (!RgX_equal(powgi(y, K), x)) return gc_long(av,0);
618 : }
619 : else
620 28 : y = pol_1(varn(x));
621 637 : if (pt)
622 : {
623 602 : if (!gequal1(a))
624 : {
625 35 : if (!b) b = gsqrtn(a, K, NULL, DEFAULTPREC);
626 35 : y = gmul(b,y);
627 : }
628 602 : if (v) y = RgX_shift_shallow(y, v);
629 602 : *pt = gc_GEN(av, y);
630 : }
631 35 : else set_avma(av);
632 637 : return 1;
633 : }
634 :
635 : long
636 1562094 : Z_ispowerall(GEN x, ulong k, GEN *pt)
637 : {
638 1562094 : long s = signe(x);
639 : ulong mask;
640 1562094 : if (!s) { if (pt) *pt = gen_0; return 1; }
641 1562094 : if (s > 0) {
642 1561765 : if (k == 2) return Z_issquareall(x, pt);
643 1425715 : if (k == 3) { mask = 1; return !!is_357_power(x, pt, &mask); }
644 228004 : if (k == 5) { mask = 2; return !!is_357_power(x, pt, &mask); }
645 60403 : if (k == 7) { mask = 4; return !!is_357_power(x, pt, &mask); }
646 54285 : return is_kth_power(x, k, pt);
647 : }
648 329 : if (!odd(k)) return 0;
649 140 : if (Z_ispowerall(absi_shallow(x), k, pt))
650 : {
651 126 : if (pt) *pt = negi(*pt);
652 126 : return 1;
653 : };
654 14 : return 0;
655 : }
656 :
657 : /* is x a K-th power mod p ? Assume p prime. */
658 : int
659 707 : Fp_ispower(GEN x, GEN K, GEN p)
660 : {
661 707 : pari_sp av = avma;
662 : GEN p_1;
663 707 : x = modii(x, p);
664 707 : if (!signe(x) || equali1(x)) return gc_bool(av,1);
665 : /* implies p > 2 */
666 224 : p_1 = subiu(p,1);
667 224 : K = gcdii(K, p_1);
668 224 : if (equali1(K)) return gc_bool(av, 1);
669 189 : if (absequaliu(K, 2)) return gc_bool(av, kronecker(x,p) > 0);
670 49 : x = Fp_pow(x, diviiexact(p_1,K), p);
671 49 : return gc_bool(av, equali1(x));
672 : }
673 :
674 : /* x unit defined modulo 2^e, e > 0, p prime */
675 : static int
676 4361 : U2_issquare(GEN x, long e)
677 : {
678 4361 : long r = signe(x)>=0?mod8(x):8-mod8(x);
679 4361 : if (e==1) return 1;
680 4347 : if (e==2) return (r&3L) == 1;
681 3549 : return r == 1;
682 : }
683 : /* x unit defined modulo p^e, e > 0, p prime */
684 : static int
685 7350 : Up_issquare(GEN x, GEN p, long e)
686 7350 : { return (absequaliu(p,2))? U2_issquare(x, e): kronecker(x,p)==1; }
687 :
688 : long
689 5558 : Zn_issquare(GEN d, GEN fn)
690 : {
691 : long j, np;
692 5558 : if (typ(d) != t_INT) pari_err_TYPE("Zn_issquare",d);
693 5558 : if (typ(fn) == t_INT) return Zn_ispower(d, fn, gen_2, NULL);
694 : /* integer factorization */
695 5439 : np = nbrows(fn);
696 9569 : for (j = 1; j <= np; ++j)
697 : {
698 9065 : GEN r, p = gcoeff(fn, j, 1);
699 9065 : long e = itos(gcoeff(fn, j, 2));
700 9065 : long v = Z_pvalrem(d,p,&r);
701 9065 : if (v < e && (odd(v) || !Up_issquare(r, p, e-v))) return 0;
702 : }
703 504 : return 1;
704 : }
705 :
706 : static long
707 1064 : Qp_ispower(GEN x, GEN K, GEN *pt)
708 : {
709 1064 : pari_sp av = avma;
710 1064 : GEN z = Qp_sqrtn(x, K, NULL);
711 1064 : if (!z) return gc_long(av,0);
712 791 : if (pt) *pt = z;
713 791 : return 1;
714 : }
715 :
716 : long
717 8558762 : ispower(GEN x, GEN K, GEN *pt)
718 : {
719 : GEN z;
720 :
721 8558762 : if (!K) return gisanypower(x, pt);
722 1558580 : if (typ(K) != t_INT) pari_err_TYPE("ispower",K);
723 1558580 : if (signe(K) <= 0) pari_err_DOMAIN("ispower","exponent","<=",gen_0,K);
724 1558580 : if (equali1(K)) { if (pt) *pt = gcopy(x); return 1; }
725 1558496 : switch(typ(x)) {
726 1485547 : case t_INT:
727 1485547 : if (lgefint(K) != 3) return 0;
728 1485491 : return Z_ispowerall(x, itou(K), pt);
729 70184 : case t_FRAC:
730 : {
731 70184 : GEN a = gel(x,1), b = gel(x,2);
732 : ulong k;
733 70184 : if (lgefint(K) != 3) return 0;
734 70177 : k = itou(K);
735 70177 : if (pt) {
736 70114 : z = cgetg(3, t_FRAC);
737 70114 : if (Z_ispowerall(a, k, &a) && Z_ispowerall(b, k, &b)) {
738 1484 : *pt = z; gel(z,1) = a; gel(z,2) = b; return 1;
739 : }
740 68630 : set_avma((pari_sp)(z + 3)); return 0;
741 : }
742 63 : return Z_ispower(a, k) && Z_ispower(b, k);
743 : }
744 609 : case t_INTMOD:
745 609 : return Zn_ispower(gel(x,2), gel(x,1), K, pt);
746 28 : case t_FFELT:
747 28 : return FF_ispower(x, K, pt);
748 :
749 1064 : case t_PADIC:
750 1064 : return Qp_ispower(x, K, pt);
751 14 : case t_POLMOD:
752 14 : return polmodispower(x, K, pt);
753 1008 : case t_POL:
754 1008 : return polispower(x, K, pt);
755 21 : case t_RFRAC:
756 21 : return rfracispower(x, K, pt);
757 7 : case t_REAL:
758 7 : if (signe(x) < 0 && !mpodd(K)) return 0;
759 : case t_COMPLEX:
760 14 : if (pt) *pt = gsqrtn(x, K, NULL, DEFAULTPREC);
761 14 : return 1;
762 :
763 7 : case t_SER:
764 7 : if (signe(x) && (!dvdsi(valser(x), K) || !ispower(gel(x,2), K, NULL)))
765 0 : return 0;
766 7 : if (pt) *pt = gsqrtn(x, K, NULL, DEFAULTPREC);
767 7 : return 1;
768 : }
769 0 : pari_err_TYPE("ispower",x);
770 : return 0; /* LCOV_EXCL_LINE */
771 : }
772 :
773 : long
774 7000182 : gisanypower(GEN x, GEN *pty)
775 : {
776 7000182 : long tx = typ(x);
777 : ulong k, h;
778 7000182 : if (tx == t_INT) return Z_isanypower(x, pty);
779 14 : if (tx == t_FRAC)
780 : {
781 14 : pari_sp av = avma;
782 14 : GEN fa, P, E, a = gel(x,1), b = gel(x,2);
783 : long i, j, p, e;
784 14 : int sw = (abscmpii(a, b) > 0);
785 :
786 14 : if (sw) swap(a, b);
787 14 : k = Z_isanypower(a, pty? &a: NULL);
788 14 : if (!k)
789 : { /* a = -1,1 or not a pure power */
790 7 : if (!is_pm1(a)) return gc_long(av,0);
791 7 : if (signe(a) < 0) b = negi(b);
792 7 : k = Z_isanypower(b, pty? &b: NULL);
793 7 : if (!k || !pty) return gc_long(av,k);
794 7 : *pty = gc_upto(av, ginv(b));
795 7 : return k;
796 : }
797 7 : fa = factoru(k);
798 7 : P = gel(fa,1);
799 7 : E = gel(fa,2); h = k;
800 14 : for (i = lg(P) - 1; i > 0; i--)
801 : {
802 7 : p = P[i];
803 7 : e = E[i];
804 21 : for (j = 0; j < e; j++)
805 14 : if (!is_kth_power(b, p, &b)) break;
806 7 : if (j < e) k /= upowuu(p, e - j);
807 : }
808 7 : if (k == 1) return gc_long(av,0);
809 7 : if (!pty) return gc_long(av,k);
810 0 : if (k != h) a = powiu(a, h/k);
811 0 : *pty = gc_GEN(av, mkfrac(a, b));
812 0 : return k;
813 : }
814 0 : pari_err_TYPE("gisanypower", x);
815 : return 0; /* LCOV_EXCL_LINE */
816 : }
817 :
818 : /* v_p(x) = e != 0 for some p; return ispower(x,,&x), updating x.
819 : * No need to optimize for 2,3,5,7 powers (done before) */
820 : static long
821 507427 : split_exponent(ulong e, GEN *x)
822 : {
823 : GEN fa, P, E;
824 507427 : long i, j, l, k = 1;
825 507427 : if (e == 1) return 1;
826 14 : fa = factoru(e);
827 14 : P = gel(fa,1);
828 14 : E = gel(fa,2); l = lg(P);
829 28 : for (i = 1; i < l; i++)
830 : {
831 14 : ulong p = P[i];
832 28 : for (j = 0; j < E[i]; j++)
833 : {
834 : GEN y;
835 14 : if (!is_kth_power(*x, p, &y)) break;
836 14 : k *= p; *x = y;
837 : }
838 : }
839 14 : return k;
840 : }
841 :
842 : /* any prime divisor of x is > 102 */
843 : static long
844 870412 : Z_isanypower_101(GEN *px)
845 : {
846 870412 : const double LOG2_103 = 6.6865; /* lower bound for log_2(103) */
847 870412 : const double LOG103 = 4.6347; /* lower bound for log(103) */
848 : forprime_t T;
849 870412 : ulong mask = 7, e2;
850 : long k, ex;
851 870412 : GEN y, x = *px;
852 :
853 870412 : k = 1;
854 873942 : while (Z_issquareall(x, &y)) { k <<= 1; x = y; }
855 870854 : while ( (ex = is_357_power(x, &y, &mask)) ) { k *= ex; x = y; }
856 870412 : e2 = (ulong)((expi(x) + 1) / LOG2_103); /* >= log_103 (x) */
857 870412 : if (u_forprime_init(&T, 11, e2))
858 : {
859 17661 : GEN logx = NULL;
860 17661 : const ulong Q = 30011; /* prime */
861 : ulong p, xmodQ;
862 17661 : double dlogx = 0;
863 : /* cut off at x^(1/p) ~ 2^30 bits which seems to be about optimum;
864 : * for large p the modular checks are no longer competitively fast */
865 17703 : while ( (ex = is_pth_power(x, &y, &T, 30)) )
866 : {
867 42 : k *= ex; x = y;
868 42 : e2 = (ulong)((expi(x) + 1) / LOG2_103);
869 42 : u_forprime_restrict(&T, e2);
870 : }
871 17661 : if (DEBUGLEVEL>4)
872 0 : err_printf("Z_isanypower: now k=%ld, x=%ld-bit\n", k, expi(x)+1);
873 17661 : xmodQ = umodiu(x, Q);
874 : /* test Q | x, just in case */
875 17661 : if (!xmodQ) { *px = x; return k * split_exponent(Z_lval(x,Q), px); }
876 : /* x^(1/p) < 2^31 */
877 17647 : p = T.p;
878 17647 : if (p <= e2)
879 : {
880 17633 : logx = logr_abs( itor(x, DEFAULTPREC) );
881 17633 : dlogx = rtodbl(logx);
882 17633 : e2 = (ulong)(dlogx / LOG103); /* >= log_103(x) */
883 : }
884 147532 : while (p && p <= e2)
885 : { /* is x a p-th power ? By computing y = round(x^(1/p)).
886 : * Check whether y^p = x, first mod Q, then exactly. */
887 129885 : pari_sp av = avma;
888 : long e;
889 129885 : GEN logy = divru(logx, p), y = grndtoi(mpexp(logy), &e);
890 129885 : ulong ymodQ = umodiu(y,Q);
891 129885 : if (e >= -10 || Fl_powu(ymodQ, p % (Q-1), Q) != xmodQ
892 129885 : || !equalii(powiu(y, p), x)) set_avma(av);
893 : else
894 : {
895 21 : k *= p; x = y; xmodQ = ymodQ; logx = logy; dlogx /= p;
896 21 : e2 = (ulong)(dlogx / LOG103); /* >= log_103(x) */
897 21 : u_forprime_restrict(&T, e2);
898 21 : continue; /* if success, retry same p */
899 : }
900 129864 : p = u_forprime_next(&T);
901 : }
902 : }
903 870398 : *px = x; return k;
904 : }
905 :
906 : static ulong tinyprimes[] = {
907 : 2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71,
908 : 73, 79, 83, 89, 97, 101, 103, 107, 109, 113, 127, 131, 137, 139, 149, 151,
909 : 157, 163, 167, 173, 179, 181, 191, 193, 197, 199
910 : };
911 :
912 : /* disregard the sign of x, caller will take care of x < 0 */
913 : static long
914 7195979 : Z_isanypower_aux(GEN x, GEN *pty)
915 : {
916 : long ex, v, i, l, k;
917 : GEN y, P, E;
918 7195979 : ulong mask, e = 0;
919 :
920 7195979 : if (abscmpii(x, gen_2) < 0) return 0; /* -1,0,1 */
921 :
922 7195249 : x = absi_shallow(x);
923 7195249 : k = l = 1;
924 7195249 : P = cgetg(26 + 1, t_VECSMALL);
925 7195249 : E = cgetg(26 + 1, t_VECSMALL);
926 : /* trial division */
927 61993288 : for(i = 0; i < 26; i++)
928 : {
929 60651366 : ulong p = tinyprimes[i];
930 : int stop;
931 60651366 : v = Z_lvalrem_stop(&x, p, &stop);
932 60651366 : if (v)
933 : {
934 8043902 : P[l] = p;
935 8043902 : E[l] = v; l++;
936 8423925 : e = ugcd(e, v); if (e == 1) goto END;
937 : }
938 55178062 : if (stop) {
939 380023 : if (is_pm1(x)) k = e;
940 380023 : goto END;
941 : }
942 : }
943 :
944 1341922 : if (e)
945 : { /* Bingo. Result divides e */
946 : long v3, v5, v7;
947 507413 : ulong e2 = e;
948 507413 : v = u_lvalrem(e2, 2, &e2);
949 507413 : if (v)
950 : {
951 377734 : for (i = 0; i < v; i++)
952 : {
953 375648 : if (!Z_issquareall(x, &y)) break;
954 2310 : k <<= 1; x = y;
955 : }
956 : }
957 507413 : mask = 0;
958 507413 : v3 = u_lvalrem(e2, 3, &e2); if (v3) mask = 1;
959 507413 : v5 = u_lvalrem(e2, 5, &e2); if (v5) mask |= 2;
960 507413 : v7 = u_lvalrem(e2, 7, &e2); if (v7) mask |= 4;
961 507497 : while ( (ex = is_357_power(x, &y, &mask)) ) {
962 84 : x = y;
963 84 : switch(ex)
964 : {
965 35 : case 3: k *= 3; if (--v3 == 0) mask &= ~1; break;
966 28 : case 5: k *= 5; if (--v5 == 0) mask &= ~2; break;
967 21 : case 7: k *= 7; if (--v7 == 0) mask &= ~4; break;
968 : }
969 : }
970 507413 : k *= split_exponent(e2, &x);
971 : }
972 : else
973 834509 : k = Z_isanypower_101(&x);
974 7195249 : END:
975 7195249 : if (pty && k != 1)
976 : {
977 70254 : if (e)
978 : { /* add missing small factors */
979 66871 : y = powuu(P[1], E[1] / k);
980 76862 : for (i = 2; i < l; i++) y = mulii(y, powuu(P[i], E[i] / k));
981 66871 : x = equali1(x)? y: mulii(x,y);
982 : }
983 70254 : *pty = x;
984 : }
985 7195249 : return k == 1? 0: k;
986 : }
987 :
988 : long
989 7195979 : Z_isanypower(GEN x, GEN *pty)
990 : {
991 7195979 : pari_sp av = avma;
992 7195979 : long k = Z_isanypower_aux(x, pty);
993 7195979 : if (!k) return gc_long(av,0);
994 70317 : if (signe(x) < 0)
995 : {
996 42 : long v = vals(k);
997 42 : if (v)
998 : {
999 28 : k >>= v;
1000 28 : if (k == 1) return gc_long(av,0);
1001 21 : if (!pty) return gc_long(av,k);
1002 14 : *pty = gc_INT(av, powiu(*pty, 1<<v));
1003 14 : togglesign(*pty); return k;
1004 : }
1005 14 : if (pty) togglesign_safe(pty);
1006 : }
1007 70289 : if (!pty) return gc_long(av, k);
1008 70233 : *pty = gc_GEN(av, *pty); return k;
1009 : }
1010 :
1011 : /* Faster than expi(n) == vali(n) or hamming(n) == 1 even for single-word
1012 : * values. If all you have is a word, you can just use n & !(n & (n-1)). */
1013 : long
1014 586618 : Z_ispow2(GEN n)
1015 : {
1016 : GEN xp;
1017 : long i, l;
1018 : ulong u;
1019 586618 : if (signe(n) != 1) return 0;
1020 530611 : xp = int_LSW(n); u = *xp; l = lgefint(n);
1021 596486 : for (i = 3; i < l; ++i)
1022 : {
1023 315550 : if (u) return 0;
1024 65875 : xp = int_nextW(xp); u = *xp;
1025 : }
1026 280936 : return !(u & (u-1));
1027 : }
1028 :
1029 : static long
1030 842200 : isprimepower_i(GEN n, GEN *pt, long flag)
1031 : {
1032 842200 : pari_sp av = avma;
1033 : long i, v;
1034 :
1035 842200 : if (typ(n) != t_INT) pari_err_TYPE("isprimepower", n);
1036 842200 : if (signe(n) <= 0) return 0;
1037 :
1038 842200 : if (lgefint(n) == 3)
1039 : {
1040 : ulong p;
1041 541260 : v = uisprimepower(n[2], &p);
1042 541260 : if (v)
1043 : {
1044 55034 : if (pt) *pt = utoipos(p);
1045 55034 : return v;
1046 : }
1047 486226 : return 0;
1048 : }
1049 1664326 : for (i = 0; i < 26; i++)
1050 : {
1051 1628423 : ulong p = tinyprimes[i];
1052 1628423 : v = Z_lvalrem(n, p, &n);
1053 1628423 : if (v)
1054 : {
1055 265037 : set_avma(av);
1056 265037 : if (!is_pm1(n)) return 0;
1057 641 : if (pt) *pt = utoipos(p);
1058 641 : return v;
1059 : }
1060 : }
1061 : /* p | n => p >= 103 */
1062 35903 : v = Z_isanypower_101(&n); /* expensive */
1063 35903 : if (!(flag? isprime(n): BPSW_psp(n))) return gc_long(av,0);
1064 5591 : if (pt) *pt = gc_GEN(av, n); else set_avma(av);
1065 5591 : return v;
1066 : }
1067 : long
1068 840147 : isprimepower(GEN n, GEN *pt) { return isprimepower_i(n,pt,1); }
1069 : long
1070 2053 : ispseudoprimepower(GEN n, GEN *pt) { return isprimepower_i(n,pt,0); }
1071 :
1072 : long
1073 546580 : uisprimepower(ulong n, ulong *pp)
1074 : { /* We must have CUTOFF^11 >= ULONG_MAX and CUTOFF^3 < ULONG_MAX.
1075 : * Tests suggest that 200-300 is the best range for 64-bit platforms. */
1076 546580 : const ulong CUTOFF = 200UL;
1077 546580 : const long TINYCUTOFF = 46; /* tinyprimes[45] = 199 */
1078 546580 : const ulong CUTOFF3 = CUTOFF*CUTOFF*CUTOFF;
1079 : #ifdef LONG_IS_64BIT
1080 : /* primes preceeding the appropriate root of ULONG_MAX. */
1081 485652 : const ulong ROOT9 = 137;
1082 485652 : const ulong ROOT8 = 251;
1083 485652 : const ulong ROOT7 = 563;
1084 485652 : const ulong ROOT5 = 7129;
1085 485652 : const ulong ROOT4 = 65521;
1086 : #else
1087 60928 : const ulong ROOT9 = 11;
1088 60928 : const ulong ROOT8 = 13;
1089 60928 : const ulong ROOT7 = 23;
1090 60928 : const ulong ROOT5 = 83;
1091 60928 : const ulong ROOT4 = 251;
1092 : #endif
1093 : ulong mask;
1094 : long v, i;
1095 : int e;
1096 546580 : if (n < 2) return 0;
1097 546531 : if (!odd(n)) {
1098 274001 : if (n & (n-1)) return 0;
1099 3658 : *pp = 2; return vals(n);
1100 : }
1101 272530 : if (n < 8) { *pp = n; return 1; } /* 3,5,7 */
1102 3654585 : for (i = 1/*skip p=2*/; i < TINYCUTOFF; i++)
1103 : {
1104 3595510 : ulong p = tinyprimes[i];
1105 3595510 : if (n % p == 0)
1106 : {
1107 212118 : v = u_lvalrem(n, p, &n);
1108 212118 : if (n == 1) { *pp = p; return v; }
1109 209977 : return 0;
1110 : }
1111 : }
1112 : /* p | n => p >= CUTOFF */
1113 :
1114 59075 : if (n < CUTOFF3)
1115 : {
1116 46354 : if (n < CUTOFF*CUTOFF || uisprime_101(n)) { *pp = n; return 1; }
1117 0 : if (uissquareall(n, &n)) { *pp = n; return 2; }
1118 0 : return 0;
1119 : }
1120 :
1121 : /* Check for squares, fourth powers, and eighth powers as appropriate. */
1122 12721 : v = 1;
1123 12721 : if (uissquareall(n, &n)) {
1124 0 : v <<= 1;
1125 0 : if (CUTOFF <= ROOT4 && uissquareall(n, &n)) {
1126 0 : v <<= 1;
1127 0 : if (CUTOFF <= ROOT8 && uissquareall(n, &n)) v <<= 1;
1128 : }
1129 : }
1130 :
1131 12721 : if (CUTOFF > ROOT5) mask = 1;
1132 : else
1133 : {
1134 12720 : const ulong CUTOFF5 = CUTOFF3*CUTOFF*CUTOFF;
1135 12720 : if (n < CUTOFF5) mask = 1; else mask = 3;
1136 12720 : if (CUTOFF <= ROOT7)
1137 : {
1138 12720 : const ulong CUTOFF7 = CUTOFF5*CUTOFF*CUTOFF;
1139 12720 : if (n >= CUTOFF7) mask = 7;
1140 : }
1141 : }
1142 :
1143 12721 : if (CUTOFF <= ROOT9 && (e = uis_357_power(n, &n, &mask))) { v *= e; mask=1; }
1144 12721 : if ((e = uis_357_power(n, &n, &mask))) v *= e;
1145 :
1146 12721 : if (uisprime_101(n)) { *pp = n; return v; }
1147 6984 : return 0;
1148 : }
1149 :
1150 : long
1151 97076 : uisnilpotent(ulong n)
1152 : {
1153 97076 : pari_sp av = avma;
1154 : long i, l;
1155 : GEN F, P, E;
1156 97076 : if (n <= 5) return 1;
1157 8533 : F = factoru(n); P = gel(F,1); E = gel(F,2);
1158 8533 : l = lg(P);
1159 8533 : if (l==2) return gc_long(av, 1);
1160 8365 : for (i=1; i < l; i++)
1161 : {
1162 8330 : ulong j, p = uel(P,i), e = uel(E,i), q = p%n;
1163 14658 : for (j = 1; j <= e; j++, q = Fl_mul(q, p, n))
1164 12292 : if (ugcd(q-1, n) > 1) return gc_long(av, 0);
1165 : }
1166 35 : return gc_long(av, 1);
1167 : }
|
