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 833184 : ulogintall(ulong B, ulong y, ulong *ptq)
29 : {
30 : ulong r, r2;
31 : long e;
32 :
33 833184 : if (y == 2)
34 : {
35 23070 : long eB = expu(B); /* 2^eB <= B < 2^(eB + 1) */
36 23070 : if (ptq) *ptq = 1UL << eB;
37 23070 : return eB;
38 : }
39 810114 : r = y, r2 = 1UL;
40 2936703 : for (e=1;; e++)
41 : { /* here, r = y^e, r2 = y^(e-1) */
42 2936703 : if (r >= B)
43 : {
44 808296 : if (r != B) { e--; r = r2; }
45 808296 : if (ptq) *ptq = r;
46 808296 : return e;
47 : }
48 2128407 : r2 = r;
49 2128407 : r = umuluu_or_0(y, r);
50 2128445 : if (!r)
51 : {
52 1856 : if (ptq) *ptq = r2;
53 1856 : 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 934108 : logintall(GEN B, GEN y, GEN *ptq)
62 : {
63 : pari_sp av;
64 934108 : long ey, e, emax, i, eB = expi(B); /* 2^eB <= B < 2^(eB + 1) */
65 : GEN q, pow2;
66 :
67 934108 : if (lgefint(B) == 3)
68 : {
69 : ulong q;
70 832034 : if (lgefint(y) > 3)
71 : {
72 0 : if (ptq) *ptq = gen_1;
73 0 : return 0;
74 : }
75 832034 : if (!ptq) return ulogintall(B[2], y[2], NULL);
76 127408 : e = ulogintall(B[2], y[2], &q);
77 127410 : *ptq = utoi(q); return e;
78 : }
79 102074 : if (equaliu(y,2))
80 : {
81 888 : if (ptq) *ptq = int2n(eB);
82 888 : return eB;
83 : }
84 101198 : av = avma;
85 101198 : ey = expi(y);
86 : /* eB/(ey+1) - 1 < e <= eB/ey */
87 101199 : emax = eB/ey;
88 101199 : if (emax <= 13) /* e small, be naive */
89 : {
90 10905 : GEN r = y, r2 = gen_1;
91 10905 : for (e=1;; e++)
92 105419 : { /* here, r = y^e, r2 = y^(e-1) */
93 116324 : long fl = cmpii(r, B);
94 116324 : if (fl >= 0)
95 : {
96 10905 : if (fl) { e--; cgiv(r); r = r2; }
97 10905 : if (ptq) *ptq = gerepileuptoint(av, r); else set_avma(av);
98 10905 : return e;
99 : }
100 105419 : 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 90294 : pow2 = new_chunk((long)log2(eB)+2);
108 90294 : gel(pow2,0) = y;
109 90294 : for (i=0, q=y;; )
110 422909 : {
111 513203 : GEN r = gel(pow2,i); /* r = y^2^i */
112 513203 : long fl = cmpii(r,B);
113 513179 : 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 513179 : if (fl > 0) { i--; break; }
120 481551 : q = r;
121 481551 : if (1L<<(i+1) > emax) break;
122 422886 : gel(pow2,++i) = sqri(q);
123 : }
124 :
125 90293 : for (e = 1L<<i;;)
126 391267 : { /* y^e = q < B < r = q * y^(2^i) */
127 481560 : pari_sp av2 = avma;
128 : long fl;
129 : GEN r;
130 481560 : if (--i < 0) break;
131 391277 : r = mulii(q, gel(pow2,i));
132 391281 : fl = cmpii(r, B);
133 391295 : if (fl > 0) set_avma(av2);
134 : else
135 : {
136 175970 : e += (1L<<i);
137 175970 : q = r;
138 175970 : if (!fl) break; /* B = r */
139 : }
140 : }
141 90290 : if (ptq) *ptq = gerepileuptoint(av, q); else set_avma(av);
142 90290 : return e;
143 : }
144 :
145 : long
146 26189 : logint0(GEN B, GEN y, GEN *ptq)
147 : {
148 26189 : const char *f = "logint";
149 26189 : if (typ(y) != t_INT) pari_err_TYPE(f,y);
150 26189 : if (cmpis(y, 2) < 0) pari_err_DOMAIN(f, "b" ,"<=", gen_1, y);
151 26189 : 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 26112 : if (signe(B) <= 0) pari_err_DOMAIN(f, "x" ,"<=", gen_0, B);
189 26112 : return logintall(B,y,ptq);
190 : }
191 :
192 : /*********************************************************************/
193 : /** INTEGRAL SQUARE ROOT **/
194 : /*********************************************************************/
195 : GEN
196 88429 : sqrtint(GEN a)
197 : {
198 88429 : 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 88373 : switch (signe(a))
225 : {
226 88352 : 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 14327240 : squaremod(ulong A)
256 : {
257 14327240 : 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 14327240 : 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 14327240 : 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 14327240 : const int squaremod11[]={1,1,0,1,1,1,0,0,0,1, 0};
267 14327240 : return (squaremod64[A & 0x3fUL]
268 6760235 : && squaremod63[A % 63UL]
269 4706373 : && squaremod65[A % 65UL]
270 21087475 : && squaremod11[A % 11UL]);
271 : }
272 :
273 : /* emulate Z_issquareall on single-word integers */
274 : long
275 10689342 : uissquareall(ulong A, ulong *sqrtA)
276 : {
277 10689342 : if (!A) { *sqrtA = 0; return 1; }
278 10689342 : if (squaremod(A))
279 : {
280 1936874 : ulong a = usqrt(A);
281 1936867 : if (a * a == A) { *sqrtA = a; return 1; }
282 : }
283 8784662 : return 0;
284 : }
285 : long
286 274831 : uissquare(ulong A)
287 : {
288 274831 : if (!A) return 1;
289 274831 : if (squaremod(A)) { ulong a = usqrt(A); if (a * a == A) return 1; }
290 265038 : return 0;
291 : }
292 :
293 : long
294 6681695 : Z_issquareall(GEN x, GEN *pt)
295 : {
296 : pari_sp av;
297 : GEN y, r;
298 :
299 6681695 : switch(signe(x))
300 : {
301 181019 : case -1: return 0;
302 1813 : case 0: if (pt) *pt=gen_0; return 1;
303 : }
304 6498863 : if (lgefint(x) == 3)
305 : {
306 3135803 : ulong u = uel(x,2), a;
307 3135803 : if (!pt) return uissquare(u);
308 2860972 : if (!uissquareall(u, &a)) return 0;
309 1476669 : *pt = utoipos(a); return 1;
310 : }
311 3363060 : if (!squaremod(umodiu(x, 64*63*65*11))) return 0;
312 2211950 : av = avma; y = sqrtremi(x, &r);
313 2211950 : if (r != gen_0) return gc_long(av,0);
314 41269 : if (pt) { *pt = y; set_avma((pari_sp)y); } else set_avma(av);
315 41269 : 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 : 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 = gerepilecopy(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 = gerepileupto(av, x); return 1;
436 : }
437 : long
438 252275 : issquareall(GEN x, GEN *pt)
439 : {
440 252275 : long tx = typ(x);
441 : GEN F;
442 : pari_sp av;
443 :
444 252275 : if (!pt) return issquare(x);
445 38093 : switch(tx)
446 : {
447 17506 : case t_INT: return Z_issquareall(x, pt);
448 1939 : case t_FRAC: av = avma;
449 1939 : F = cgetg(3, t_FRAC);
450 1939 : if ( !Z_issquareall(gel(x,1), &gel(F,1))
451 1939 : || !Z_issquareall(gel(x,2), &gel(F,2))) return gc_long(av,0);
452 1876 : *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_ispower(gel(x,2), gel(x,1), gen_2, NULL);
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 = gel(x,4); if (!signe(a)) return 1;
500 133 : if (valp(x)&1) return 0;
501 119 : p = gel(x,2);
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 0 : GEN gissquare(GEN x) { return issquare(x)? gen_1: gen_0; }
527 0 : GEN gissquareall(GEN x, GEN *pt) { return issquareall(x,pt)? gen_1: gen_0; }
528 :
529 : long
530 1386 : ispolygonal(GEN x, GEN S, GEN *N)
531 : {
532 1386 : pari_sp av = avma;
533 : GEN D, d, n;
534 1386 : if (typ(x) != t_INT) pari_err_TYPE("ispolygonal", x);
535 1386 : if (typ(S) != t_INT) pari_err_TYPE("ispolygonal", S);
536 1386 : if (abscmpiu(S,3) < 0) pari_err_DOMAIN("ispolygonal","s","<", utoipos(3),S);
537 1386 : if (signe(x) < 0) return 0;
538 1386 : if (signe(x) == 0) { if (N) *N = gen_0; return 1; }
539 1260 : if (is_pm1(x)) { if (N) *N = gen_1; return 1; }
540 : /* n = (sqrt( (8s - 16) x + (s-4)^2 ) + s - 4) / 2(s - 2) */
541 1134 : if (abscmpiu(S, 1<<16) < 0) /* common case ! */
542 : {
543 441 : ulong s = S[2], r;
544 441 : if (s == 4) return Z_issquareall(x, N);
545 378 : if (s == 3)
546 0 : D = addiu(shifti(x, 3), 1);
547 : else
548 378 : D = addiu(mului(8*s - 16, x), (s-4)*(s-4));
549 378 : if (!Z_issquareall(D, &d)) return gc_long(av,0);
550 378 : if (s == 3)
551 0 : d = subiu(d, 1);
552 : else
553 378 : d = addiu(d, s - 4);
554 378 : n = absdiviu_rem(d, 2*s - 4, &r);
555 378 : if (r) return gc_long(av,0);
556 : }
557 : else
558 : {
559 693 : GEN r, S_2 = subiu(S,2), S_4 = subiu(S,4);
560 693 : D = addii(mulii(shifti(S_2,3), x), sqri(S_4));
561 693 : if (!Z_issquareall(D, &d)) return gc_long(av,0);
562 693 : d = addii(d, S_4);
563 693 : n = dvmdii(shifti(d,-1), S_2, &r);
564 693 : if (r != gen_0) return gc_long(av,0);
565 : }
566 1071 : if (N) *N = gerepileuptoint(av, n); else set_avma(av);
567 1071 : return 1;
568 : }
569 :
570 : /*********************************************************************/
571 : /** PERFECT POWER **/
572 : /*********************************************************************/
573 : static long
574 987 : polispower(GEN x, GEN K, GEN *pt)
575 : {
576 : pari_sp av;
577 987 : long v, d, k = itos(K);
578 : GEN y, a, b;
579 987 : GEN T = NULL, p = NULL;
580 :
581 987 : if (!signe(x))
582 : {
583 7 : if (pt) *pt = gcopy(x);
584 7 : return 1;
585 : }
586 980 : d = degpol(x);
587 980 : if (d % k) return 0; /* degree not multiple of k */
588 973 : av = avma;
589 973 : if (RgX_is_FpXQX(x, &T, &p) && p)
590 : { /* over Fq */
591 336 : if (T && typ(T) == t_FFELT)
592 : {
593 126 : if (!FFX_ispower(x, k, T, pt)) return gc_long(av,0);
594 105 : return 1;
595 : }
596 210 : x = RgX_to_FqX(x,T,p);
597 210 : if (!FqX_ispower(x, k, T,p, pt)) return gc_long(av,0);
598 175 : if (pt) *pt = gerepileupto(av, FqX_to_mod(*pt, T, p));
599 175 : return 1;
600 : }
601 637 : v = RgX_valrem(x, &x);
602 637 : if (v % k) return 0;
603 630 : v /= k;
604 630 : a = gel(x,2); b = NULL;
605 630 : if (!ispower(a, K, &b)) return gc_long(av,0);
606 616 : if (d)
607 : {
608 588 : GEN p = characteristic(x);
609 588 : a = leading_coeff(x);
610 588 : if (!ispower(a, K, &b)) return gc_long(av,0);
611 588 : x = RgX_normalize(x);
612 588 : if (signe(p) && cmpii(p,K) <= 0)
613 0 : pari_err_IMPL("ispower(general t_POL) in small characteristic");
614 588 : y = gtrunc(gsqrtn(RgX_to_ser(x,lg(x)), K, NULL, 0));
615 588 : if (!RgX_equal(powgi(y, K), x)) return gc_long(av,0);
616 : }
617 : else
618 28 : y = pol_1(varn(x));
619 616 : if (pt)
620 : {
621 581 : if (!gequal1(a))
622 : {
623 35 : if (!b) b = gsqrtn(a, K, NULL, DEFAULTPREC);
624 35 : y = gmul(b,y);
625 : }
626 581 : if (v) y = RgX_shift_shallow(y, v);
627 581 : *pt = gerepilecopy(av, y);
628 : }
629 35 : else set_avma(av);
630 616 : return 1;
631 : }
632 :
633 : long
634 1559328 : Z_ispowerall(GEN x, ulong k, GEN *pt)
635 : {
636 1559328 : long s = signe(x);
637 : ulong mask;
638 1559328 : if (!s) { if (pt) *pt = gen_0; return 1; }
639 1559328 : if (s > 0) {
640 1558999 : if (k == 2) return Z_issquareall(x, pt);
641 1425582 : if (k == 3) { mask = 1; return !!is_357_power(x, pt, &mask); }
642 227878 : if (k == 5) { mask = 2; return !!is_357_power(x, pt, &mask); }
643 60305 : if (k == 7) { mask = 4; return !!is_357_power(x, pt, &mask); }
644 54187 : return is_kth_power(x, k, pt);
645 : }
646 329 : if (!odd(k)) return 0;
647 140 : if (Z_ispowerall(absi_shallow(x), k, pt))
648 : {
649 126 : if (pt) *pt = negi(*pt);
650 126 : return 1;
651 : };
652 14 : return 0;
653 : }
654 :
655 : /* is x a K-th power mod p ? Assume p prime. */
656 : int
657 567 : Fp_ispower(GEN x, GEN K, GEN p)
658 : {
659 567 : pari_sp av = avma;
660 : GEN p_1;
661 567 : x = modii(x, p);
662 567 : if (!signe(x) || equali1(x)) return gc_bool(av,1);
663 : /* implies p > 2 */
664 126 : p_1 = subiu(p,1);
665 126 : K = gcdii(K, p_1);
666 126 : if (absequaliu(K, 2)) return gc_bool(av, kronecker(x,p) > 0);
667 49 : x = Fp_pow(x, diviiexact(p_1,K), p);
668 49 : return gc_bool(av, equali1(x));
669 : }
670 :
671 : /* x unit defined modulo 2^e, e > 0, p prime */
672 : static int
673 2541 : U2_issquare(GEN x, long e)
674 : {
675 2541 : long r = signe(x)>=0?mod8(x):8-mod8(x);
676 2541 : if (e==1) return 1;
677 2541 : if (e==2) return (r&3L) == 1;
678 2009 : return r == 1;
679 : }
680 : /* x unit defined modulo p^e, e > 0, p prime */
681 : static int
682 5229 : Up_issquare(GEN x, GEN p, long e)
683 5229 : { return (absequaliu(p,2))? U2_issquare(x, e): kronecker(x,p)==1; }
684 :
685 : long
686 2793 : Zn_issquare(GEN d, GEN fn)
687 : {
688 : long j, np;
689 2793 : if (typ(d) != t_INT) pari_err_TYPE("Zn_issquare",d);
690 2793 : if (typ(fn) == t_INT) return Zn_ispower(d, fn, gen_2, NULL);
691 : /* integer factorization */
692 2793 : np = nbrows(fn);
693 6006 : for (j = 1; j <= np; ++j)
694 : {
695 5586 : GEN r, p = gcoeff(fn, j, 1);
696 5586 : long e = itos(gcoeff(fn, j, 2));
697 5586 : long v = Z_pvalrem(d,p,&r);
698 5586 : if (v < e && (odd(v) || !Up_issquare(r, p, e-v))) return 0;
699 : }
700 420 : return 1;
701 : }
702 :
703 : static long
704 1113 : Qp_ispower(GEN x, GEN K, GEN *pt)
705 : {
706 1113 : pari_sp av = avma;
707 1113 : GEN z = Qp_sqrtn(x, K, NULL);
708 1113 : if (!z) return gc_long(av,0);
709 819 : if (pt) *pt = z;
710 819 : return 1;
711 : }
712 :
713 : long
714 8556102 : ispower(GEN x, GEN K, GEN *pt)
715 : {
716 : GEN z;
717 :
718 8556102 : if (!K) return gisanypower(x, pt);
719 1555920 : if (typ(K) != t_INT) pari_err_TYPE("ispower",K);
720 1555920 : if (signe(K) <= 0) pari_err_DOMAIN("ispower","exponent","<=",gen_0,K);
721 1555920 : if (equali1(K)) { if (pt) *pt = gcopy(x); return 1; }
722 1555857 : switch(typ(x)) {
723 1482880 : case t_INT:
724 1482880 : if (lgefint(K) != 3) return 0;
725 1482824 : return Z_ispowerall(x, itou(K), pt);
726 70184 : case t_FRAC:
727 : {
728 70184 : GEN a = gel(x,1), b = gel(x,2);
729 : ulong k;
730 70184 : if (lgefint(K) != 3) return 0;
731 70177 : k = itou(K);
732 70177 : if (pt) {
733 70114 : z = cgetg(3, t_FRAC);
734 70114 : if (Z_ispowerall(a, k, &a) && Z_ispowerall(b, k, &b)) {
735 1484 : *pt = z; gel(z,1) = a; gel(z,2) = b; return 1;
736 : }
737 68630 : set_avma((pari_sp)(z + 3)); return 0;
738 : }
739 63 : return Z_ispower(a, k) && Z_ispower(b, k);
740 : }
741 609 : case t_INTMOD:
742 609 : return Zn_ispower(gel(x,2), gel(x,1), K, pt);
743 28 : case t_FFELT:
744 28 : return FF_ispower(x, K, pt);
745 :
746 1113 : case t_PADIC:
747 1113 : return Qp_ispower(x, K, pt);
748 14 : case t_POLMOD:
749 14 : return polmodispower(x, K, pt);
750 987 : case t_POL:
751 987 : return polispower(x, K, pt);
752 21 : case t_RFRAC:
753 21 : return rfracispower(x, K, pt);
754 7 : case t_REAL:
755 7 : if (signe(x) < 0 && !mpodd(K)) return 0;
756 : case t_COMPLEX:
757 14 : if (pt) *pt = gsqrtn(x, K, NULL, DEFAULTPREC);
758 14 : return 1;
759 :
760 7 : case t_SER:
761 7 : if (signe(x) && (!dvdsi(valser(x), K) || !ispower(gel(x,2), K, NULL)))
762 0 : return 0;
763 7 : if (pt) *pt = gsqrtn(x, K, NULL, DEFAULTPREC);
764 7 : return 1;
765 : }
766 0 : pari_err_TYPE("ispower",x);
767 : return 0; /* LCOV_EXCL_LINE */
768 : }
769 :
770 : long
771 7000182 : gisanypower(GEN x, GEN *pty)
772 : {
773 7000182 : long tx = typ(x);
774 : ulong k, h;
775 7000182 : if (tx == t_INT) return Z_isanypower(x, pty);
776 14 : if (tx == t_FRAC)
777 : {
778 14 : pari_sp av = avma;
779 14 : GEN fa, P, E, a = gel(x,1), b = gel(x,2);
780 : long i, j, p, e;
781 14 : int sw = (abscmpii(a, b) > 0);
782 :
783 14 : if (sw) swap(a, b);
784 14 : k = Z_isanypower(a, pty? &a: NULL);
785 14 : if (!k)
786 : { /* a = -1,1 or not a pure power */
787 7 : if (!is_pm1(a)) return gc_long(av,0);
788 7 : if (signe(a) < 0) b = negi(b);
789 7 : k = Z_isanypower(b, pty? &b: NULL);
790 7 : if (!k || !pty) return gc_long(av,k);
791 7 : *pty = gerepileupto(av, ginv(b));
792 7 : return k;
793 : }
794 7 : fa = factoru(k);
795 7 : P = gel(fa,1);
796 7 : E = gel(fa,2); h = k;
797 14 : for (i = lg(P) - 1; i > 0; i--)
798 : {
799 7 : p = P[i];
800 7 : e = E[i];
801 21 : for (j = 0; j < e; j++)
802 14 : if (!is_kth_power(b, p, &b)) break;
803 7 : if (j < e) k /= upowuu(p, e - j);
804 : }
805 7 : if (k == 1) return gc_long(av,0);
806 7 : if (!pty) return gc_long(av,k);
807 0 : if (k != h) a = powiu(a, h/k);
808 0 : *pty = gerepilecopy(av, mkfrac(a, b));
809 0 : return k;
810 : }
811 0 : pari_err_TYPE("gisanypower", x);
812 : return 0; /* LCOV_EXCL_LINE */
813 : }
814 :
815 : /* v_p(x) = e != 0 for some p; return ispower(x,,&x), updating x.
816 : * No need to optimize for 2,3,5,7 powers (done before) */
817 : static long
818 507443 : split_exponent(ulong e, GEN *x)
819 : {
820 : GEN fa, P, E;
821 507443 : long i, j, l, k = 1;
822 507443 : if (e == 1) return 1;
823 14 : fa = factoru(e);
824 14 : P = gel(fa,1);
825 14 : E = gel(fa,2); l = lg(P);
826 28 : for (i = 1; i < l; i++)
827 : {
828 14 : ulong p = P[i];
829 28 : for (j = 0; j < E[i]; j++)
830 : {
831 : GEN y;
832 14 : if (!is_kth_power(*x, p, &y)) break;
833 14 : k *= p; *x = y;
834 : }
835 : }
836 14 : return k;
837 : }
838 :
839 : /* any prime divisor of x is > 102 */
840 : static long
841 870335 : Z_isanypower_101(GEN *px)
842 : {
843 870335 : const double LOG2_103 = 6.6865; /* lower bound for log_2(103) */
844 870335 : const double LOG103 = 4.6347; /* lower bound for log(103) */
845 : forprime_t T;
846 870335 : ulong mask = 7, e2;
847 : long k, ex;
848 870335 : GEN y, x = *px;
849 :
850 870335 : k = 1;
851 873837 : while (Z_issquareall(x, &y)) { k <<= 1; x = y; }
852 870770 : while ( (ex = is_357_power(x, &y, &mask)) ) { k *= ex; x = y; }
853 870335 : e2 = (ulong)((expi(x) + 1) / LOG2_103); /* >= log_103 (x) */
854 870335 : if (u_forprime_init(&T, 11, e2))
855 : {
856 17647 : GEN logx = NULL;
857 17647 : const ulong Q = 30011; /* prime */
858 : ulong p, xmodQ;
859 17647 : double dlogx = 0;
860 : /* cut off at x^(1/p) ~ 2^30 bits which seems to be about optimum;
861 : * for large p the modular checks are no longer competitively fast */
862 17689 : while ( (ex = is_pth_power(x, &y, &T, 30)) )
863 : {
864 42 : k *= ex; x = y;
865 42 : e2 = (ulong)((expi(x) + 1) / LOG2_103);
866 42 : u_forprime_restrict(&T, e2);
867 : }
868 17647 : if (DEBUGLEVEL>4)
869 0 : err_printf("Z_isanypower: now k=%ld, x=%ld-bit\n", k, expi(x)+1);
870 17647 : xmodQ = umodiu(x, Q);
871 : /* test Q | x, just in case */
872 17647 : if (!xmodQ) { *px = x; return k * split_exponent(Z_lval(x,Q), px); }
873 : /* x^(1/p) < 2^31 */
874 17633 : p = T.p;
875 17633 : if (p <= e2)
876 : {
877 17619 : logx = logr_abs( itor(x, DEFAULTPREC) );
878 17619 : dlogx = rtodbl(logx);
879 17619 : e2 = (ulong)(dlogx / LOG103); /* >= log_103(x) */
880 : }
881 147329 : while (p && p <= e2)
882 : { /* is x a p-th power ? By computing y = round(x^(1/p)).
883 : * Check whether y^p = x, first mod Q, then exactly. */
884 129696 : pari_sp av = avma;
885 : long e;
886 129696 : GEN logy = divru(logx, p), y = grndtoi(mpexp(logy), &e);
887 129696 : ulong ymodQ = umodiu(y,Q);
888 129696 : if (e >= -10 || Fl_powu(ymodQ, p % (Q-1), Q) != xmodQ
889 129696 : || !equalii(powiu(y, p), x)) set_avma(av);
890 : else
891 : {
892 21 : k *= p; x = y; xmodQ = ymodQ; logx = logy; dlogx /= p;
893 21 : e2 = (ulong)(dlogx / LOG103); /* >= log_103(x) */
894 21 : u_forprime_restrict(&T, e2);
895 21 : continue; /* if success, retry same p */
896 : }
897 129675 : p = u_forprime_next(&T);
898 : }
899 : }
900 870321 : *px = x; return k;
901 : }
902 :
903 : static ulong tinyprimes[] = {
904 : 2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71,
905 : 73, 79, 83, 89, 97, 101, 103, 107, 109, 113, 127, 131, 137, 139, 149, 151,
906 : 157, 163, 167, 173, 179, 181, 191, 193, 197, 199
907 : };
908 :
909 : /* disregard the sign of x, caller will take care of x < 0 */
910 : static long
911 7195980 : Z_isanypower_aux(GEN x, GEN *pty)
912 : {
913 : long ex, v, i, l, k;
914 : GEN y, P, E;
915 7195980 : ulong mask, e = 0;
916 :
917 7195980 : if (abscmpii(x, gen_2) < 0) return 0; /* -1,0,1 */
918 :
919 7195251 : x = absi_shallow(x);
920 7195251 : k = l = 1;
921 7195251 : P = cgetg(26 + 1, t_VECSMALL);
922 7195251 : E = cgetg(26 + 1, t_VECSMALL);
923 : /* trial division */
924 61992459 : for(i = 0; i < 26; i++)
925 : {
926 60650577 : ulong p = tinyprimes[i];
927 : int stop;
928 60650577 : v = Z_lvalrem_stop(&x, p, &stop);
929 60650576 : if (v)
930 : {
931 8043941 : P[l] = p;
932 8043941 : E[l] = v; l++;
933 8423997 : e = ugcd(e, v); if (e == 1) goto END;
934 : }
935 55177263 : if (stop) {
936 380056 : if (is_pm1(x)) k = e;
937 380056 : goto END;
938 : }
939 : }
940 :
941 1341882 : if (e)
942 : { /* Bingo. Result divides e */
943 : long v3, v5, v7;
944 507429 : ulong e2 = e;
945 507429 : v = u_lvalrem(e2, 2, &e2);
946 507429 : if (v)
947 : {
948 377780 : for (i = 0; i < v; i++)
949 : {
950 375659 : if (!Z_issquareall(x, &y)) break;
951 2333 : k <<= 1; x = y;
952 : }
953 : }
954 507429 : mask = 0;
955 507429 : v3 = u_lvalrem(e2, 3, &e2); if (v3) mask = 1;
956 507429 : v5 = u_lvalrem(e2, 5, &e2); if (v5) mask |= 2;
957 507429 : v7 = u_lvalrem(e2, 7, &e2); if (v7) mask |= 4;
958 507506 : while ( (ex = is_357_power(x, &y, &mask)) ) {
959 77 : x = y;
960 77 : switch(ex)
961 : {
962 28 : case 3: k *= 3; if (--v3 == 0) mask &= ~1; break;
963 28 : case 5: k *= 5; if (--v5 == 0) mask &= ~2; break;
964 21 : case 7: k *= 7; if (--v7 == 0) mask &= ~4; break;
965 : }
966 : }
967 507429 : k *= split_exponent(e2, &x);
968 : }
969 : else
970 834453 : k = Z_isanypower_101(&x);
971 7195251 : END:
972 7195251 : if (pty && k != 1)
973 : {
974 70224 : if (e)
975 : { /* add missing small factors */
976 66869 : y = powuu(P[1], E[1] / k);
977 76863 : for (i = 2; i < l; i++) y = mulii(y, powuu(P[i], E[i] / k));
978 66869 : x = equali1(x)? y: mulii(x,y);
979 : }
980 70224 : *pty = x;
981 : }
982 7195251 : return k == 1? 0: k;
983 : }
984 :
985 : long
986 7195980 : Z_isanypower(GEN x, GEN *pty)
987 : {
988 7195980 : pari_sp av = avma;
989 7195980 : long k = Z_isanypower_aux(x, pty);
990 7195980 : if (!k) return gc_long(av,0);
991 70287 : if (signe(x) < 0)
992 : {
993 42 : long v = vals(k);
994 42 : if (v)
995 : {
996 28 : k >>= v;
997 28 : if (k == 1) return gc_long(av,0);
998 21 : if (!pty) return gc_long(av,k);
999 14 : *pty = gerepileuptoint(av, powiu(*pty, 1<<v));
1000 14 : togglesign(*pty); return k;
1001 : }
1002 14 : if (pty) togglesign_safe(pty);
1003 : }
1004 70259 : if (!pty) return gc_long(av, k);
1005 70203 : *pty = gerepilecopy(av, *pty); return k;
1006 : }
1007 :
1008 : /* Faster than expi(n) == vali(n) or hamming(n) == 1 even for single-word
1009 : * values. If all you have is a word, you can just use n & !(n & (n-1)). */
1010 : long
1011 526511 : Z_ispow2(GEN n)
1012 : {
1013 : GEN xp;
1014 : long i, l;
1015 : ulong u;
1016 526511 : if (signe(n) != 1) return 0;
1017 470504 : xp = int_LSW(n); u = *xp; l = lgefint(n);
1018 531385 : for (i = 3; i < l; ++i)
1019 : {
1020 250463 : if (u) return 0;
1021 60881 : xp = int_nextW(xp); u = *xp;
1022 : }
1023 280922 : return !(u & (u-1));
1024 : }
1025 :
1026 : static long
1027 842128 : isprimepower_i(GEN n, GEN *pt, long flag)
1028 : {
1029 842128 : pari_sp av = avma;
1030 : long i, v;
1031 :
1032 842128 : if (typ(n) != t_INT) pari_err_TYPE("isprimepower", n);
1033 842128 : if (signe(n) <= 0) return 0;
1034 :
1035 842128 : if (lgefint(n) == 3)
1036 : {
1037 : ulong p;
1038 541232 : v = uisprimepower(n[2], &p);
1039 541232 : if (v)
1040 : {
1041 55006 : if (pt) *pt = utoipos(p);
1042 55006 : return v;
1043 : }
1044 486226 : return 0;
1045 : }
1046 1663673 : for (i = 0; i < 26; i++)
1047 : {
1048 1627791 : ulong p = tinyprimes[i];
1049 1627791 : v = Z_lvalrem(n, p, &n);
1050 1627790 : if (v)
1051 : {
1052 265013 : set_avma(av);
1053 265013 : if (!is_pm1(n)) return 0;
1054 617 : if (pt) *pt = utoipos(p);
1055 617 : return v;
1056 : }
1057 : }
1058 : /* p | n => p >= 103 */
1059 35882 : v = Z_isanypower_101(&n); /* expensive */
1060 35882 : if (!(flag? isprime(n): BPSW_psp(n))) return gc_long(av,0);
1061 5570 : if (pt) *pt = gerepilecopy(av, n); else set_avma(av);
1062 5570 : return v;
1063 : }
1064 : long
1065 840147 : isprimepower(GEN n, GEN *pt) { return isprimepower_i(n,pt,1); }
1066 : long
1067 1981 : ispseudoprimepower(GEN n, GEN *pt) { return isprimepower_i(n,pt,0); }
1068 :
1069 : long
1070 642976 : uisprimepower(ulong n, ulong *pp)
1071 : { /* We must have CUTOFF^11 >= ULONG_MAX and CUTOFF^3 < ULONG_MAX.
1072 : * Tests suggest that 200-300 is the best range for 64-bit platforms. */
1073 642976 : const ulong CUTOFF = 200UL;
1074 642976 : const long TINYCUTOFF = 46; /* tinyprimes[45] = 199 */
1075 642976 : const ulong CUTOFF3 = CUTOFF*CUTOFF*CUTOFF;
1076 : #ifdef LONG_IS_64BIT
1077 : /* primes preceeding the appropriate root of ULONG_MAX. */
1078 568277 : const ulong ROOT9 = 137;
1079 568277 : const ulong ROOT8 = 251;
1080 568277 : const ulong ROOT7 = 563;
1081 568277 : const ulong ROOT5 = 7129;
1082 568277 : const ulong ROOT4 = 65521;
1083 : #else
1084 74699 : const ulong ROOT9 = 11;
1085 74699 : const ulong ROOT8 = 13;
1086 74699 : const ulong ROOT7 = 23;
1087 74699 : const ulong ROOT5 = 83;
1088 74699 : const ulong ROOT4 = 251;
1089 : #endif
1090 : ulong mask;
1091 : long v, i;
1092 : int e;
1093 642976 : if (n < 2) return 0;
1094 642934 : if (!odd(n)) {
1095 341851 : if (n & (n-1)) return 0;
1096 65544 : *pp = 2; return vals(n);
1097 : }
1098 301083 : if (n < 8) { *pp = n; return 1; } /* 3,5,7 */
1099 3655110 : for (i = 1/*skip p=2*/; i < TINYCUTOFF; i++)
1100 : {
1101 3596035 : ulong p = tinyprimes[i];
1102 3596035 : if (n % p == 0)
1103 : {
1104 212447 : v = u_lvalrem(n, p, &n);
1105 212447 : if (n == 1) { *pp = p; return v; }
1106 210012 : return 0;
1107 : }
1108 : }
1109 : /* p | n => p >= CUTOFF */
1110 :
1111 59075 : if (n < CUTOFF3)
1112 : {
1113 46354 : if (n < CUTOFF*CUTOFF || uisprime_101(n)) { *pp = n; return 1; }
1114 0 : if (uissquareall(n, &n)) { *pp = n; return 2; }
1115 0 : return 0;
1116 : }
1117 :
1118 : /* Check for squares, fourth powers, and eighth powers as appropriate. */
1119 12721 : v = 1;
1120 12721 : if (uissquareall(n, &n)) {
1121 0 : v <<= 1;
1122 0 : if (CUTOFF <= ROOT4 && uissquareall(n, &n)) {
1123 0 : v <<= 1;
1124 0 : if (CUTOFF <= ROOT8 && uissquareall(n, &n)) v <<= 1;
1125 : }
1126 : }
1127 :
1128 12721 : if (CUTOFF > ROOT5) mask = 1;
1129 : else
1130 : {
1131 12720 : const ulong CUTOFF5 = CUTOFF3*CUTOFF*CUTOFF;
1132 12720 : if (n < CUTOFF5) mask = 1; else mask = 3;
1133 12720 : if (CUTOFF <= ROOT7)
1134 : {
1135 12720 : const ulong CUTOFF7 = CUTOFF5*CUTOFF*CUTOFF;
1136 12720 : if (n >= CUTOFF7) mask = 7;
1137 : }
1138 : }
1139 :
1140 12721 : if (CUTOFF <= ROOT9 && (e = uis_357_power(n, &n, &mask))) { v *= e; mask=1; }
1141 12721 : if ((e = uis_357_power(n, &n, &mask))) v *= e;
1142 :
1143 12721 : if (uisprime_101(n)) { *pp = n; return v; }
1144 6984 : return 0;
1145 : }
1146 :
|
