Code coverage tests

This page documents the degree to which the PARI/GP source code is tested by our public test suite, distributed with the source distribution in directory src/test/. This is measured by the gcov utility; we then process gcov output using the lcov frond-end.

We test a few variants depending on Configure flags on the pari.math.u-bordeaux.fr machine (x86_64 architecture), and agregate them in the final report:

The target is to exceed 90% coverage for all mathematical modules (given that branches depending on DEBUGLEVEL or DEBUGMEM are not covered). This script is run to produce the results below.

LCOV - code coverage report
Current view: top level - basemath - ispower.c (source / functions) Hit Total Coverage
Test: PARI/GP v2.16.2 lcov report (development 29115-f22e516b23) Lines: 686 736 93.2 %
Date: 2024-03-28 08:06:56 Functions: 37 39 94.9 %
Legend: Lines: hit not hit

          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      820453 : ulogintall(ulong B, ulong y, ulong *ptq)
      29             : {
      30             :   ulong r, r2;
      31             :   long e;
      32             : 
      33      820453 :   if (y == 2)
      34             :   {
      35       22033 :     long eB = expu(B); /* 2^eB <= B < 2^(eB + 1) */
      36       22033 :     if (ptq) *ptq = 1UL << eB;
      37       22033 :     return eB;
      38             :   }
      39      798420 :   r = y, r2 = 1UL;
      40     2753679 :   for (e=1;; e++)
      41             :   { /* here, r = y^e, r2 = y^(e-1) */
      42     2753679 :     if (r >= B)
      43             :     {
      44      796644 :       if (r != B) { e--; r = r2; }
      45      796644 :       if (ptq) *ptq = r;
      46      796644 :       return e;
      47             :     }
      48     1957035 :     r2 = r;
      49     1957035 :     r = umuluu_or_0(y, r);
      50     1957078 :     if (!r)
      51             :     {
      52        1819 :       if (ptq) *ptq = r2;
      53        1819 :       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      918234 : logintall(GEN B, GEN y, GEN *ptq)
      62             : {
      63             :   pari_sp av;
      64      918234 :   long ey, e, emax, i, eB = expi(B); /* 2^eB <= B < 2^(eB + 1) */
      65             :   GEN q, pow2;
      66             : 
      67      918237 :   if (lgefint(B) == 3)
      68             :   {
      69             :     ulong q;
      70      819303 :     if (lgefint(y) > 3)
      71             :     {
      72           0 :       if (ptq) *ptq = gen_1;
      73           0 :       return 0;
      74             :     }
      75      819303 :     if (!ptq) return ulogintall(B[2], y[2], NULL);
      76      121678 :     e = ulogintall(B[2], y[2], &q);
      77      121681 :     *ptq = utoi(q); return e;
      78             :   }
      79       98934 :   if (equaliu(y,2))
      80             :   {
      81         790 :     if (ptq) *ptq = int2n(eB);
      82         790 :     return eB;
      83             :   }
      84       98161 :   av = avma;
      85       98161 :   ey = expi(y);
      86             :   /* eB/(ey+1) - 1 < e <= eB/ey */
      87       98161 :   emax = eB/ey;
      88       98161 :   if (emax <= 13) /* e small, be naive */
      89             :   {
      90       10832 :     GEN r = y, r2 = gen_1;
      91       10832 :     for (e=1;; e++)
      92      104705 :     { /* here, r = y^e, r2 = y^(e-1) */
      93      115537 :       long fl = cmpii(r, B);
      94      115537 :       if (fl >= 0)
      95             :       {
      96       10832 :         if (fl) { e--; cgiv(r); r = r2; }
      97       10832 :         if (ptq) *ptq = gerepileuptoint(av, r); else set_avma(av);
      98       10832 :         return e;
      99             :       }
     100      104705 :       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       87329 :   pow2 = new_chunk((long)log2(eB)+2);
     108       87330 :   gel(pow2,0) = y;
     109       87330 :   for (i=0, q=y;; )
     110      408784 :   {
     111      496114 :     GEN r = gel(pow2,i); /* r = y^2^i */
     112      496114 :     long fl = cmpii(r,B);
     113      496090 :     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      496090 :     if (fl > 0) { i--; break; }
     120      465155 :     q = r;
     121      465155 :     if (1L<<(i+1) > emax) break;
     122      408762 :     gel(pow2,++i) = sqri(q);
     123             :   }
     124             : 
     125       87328 :   for (e = 1L<<i;;)
     126      377850 :   { /* y^e = q < B < r = q * y^(2^i) */
     127      465178 :     pari_sp av2 = avma;
     128             :     long fl;
     129             :     GEN r;
     130      465178 :     if (--i < 0) break;
     131      377860 :     r = mulii(q, gel(pow2,i));
     132      377867 :     fl = cmpii(r, B);
     133      377897 :     if (fl > 0) set_avma(av2);
     134             :     else
     135             :     {
     136      170693 :       e += (1L<<i);
     137      170693 :       q = r;
     138      170693 :       if (!fl) break; /* B = r */
     139             :     }
     140             :   }
     141       87325 :   if (ptq) *ptq = gerepileuptoint(av, q); else set_avma(av);
     142       87325 :   return e;
     143             : }
     144             : 
     145             : long
     146       27153 : logint0(GEN B, GEN y, GEN *ptq)
     147             : {
     148       27153 :   const char *f = "logint";
     149       27153 :   if (typ(y) != t_INT) pari_err_TYPE(f,y);
     150       27153 :   if (cmpis(y, 2) < 0) pari_err_DOMAIN(f, "b" ,"<=", gen_1, y);
     151       27153 :   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       27076 :   if (signe(B) <= 0) pari_err_DOMAIN(f, "x" ,"<=", gen_0, B);
     189       27076 :   return logintall(B,y,ptq);
     190             : }
     191             : 
     192             : /*********************************************************************/
     193             : /**                     INTEGRAL SQUARE ROOT                        **/
     194             : /*********************************************************************/
     195             : GEN
     196       87589 : sqrtint(GEN a)
     197             : {
     198       87589 :   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       87533 :   switch (signe(a))
     225             :   {
     226       87512 :     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    20411911 : squaremod(ulong A)
     256             : {
     257    20411911 :   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    20411911 :   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    20411911 :   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    20411911 :   const int squaremod11[]={1,1,0,1,1,1,0,0,0,1, 0};
     267    20411911 :   return (squaremod64[A & 0x3fUL]
     268     8017170 :     && squaremod63[A % 63UL]
     269     4872529 :     && squaremod65[A % 65UL]
     270    28429081 :     && squaremod11[A % 11UL]);
     271             : }
     272             : 
     273             : /* emulate Z_issquareall on single-word integers */
     274             : long
     275    16775836 : uissquareall(ulong A, ulong *sqrtA)
     276             : {
     277    16775836 :   if (!A) { *sqrtA = 0; return 1; }
     278    16775836 :   if (squaremod(A))
     279             :   {
     280     1726125 :     ulong a = usqrt(A);
     281     1726114 :     if (a * a == A) { *sqrtA = a; return 1; }
     282             :   }
     283    15166942 :   return 0;
     284             : }
     285             : long
     286      274286 : uissquare(ulong A)
     287             : {
     288      274286 :   if (!A) return 1;
     289      274286 :   if (squaremod(A)) { ulong a = usqrt(A); if (a * a == A) return 1; }
     290      264562 :   return 0;
     291             : }
     292             : 
     293             : long
     294     6376080 : Z_issquareall(GEN x, GEN *pt)
     295             : {
     296             :   pari_sp av;
     297             :   GEN y, r;
     298             : 
     299     6376080 :   switch(signe(x))
     300             :   {
     301      180543 :     case -1: return 0;
     302        1806 :     case 0: if (pt) *pt=gen_0; return 1;
     303             :   }
     304     6193731 :   if (lgefint(x) == 3)
     305             :   {
     306     2831945 :     ulong u = uel(x,2), a;
     307     2831945 :     if (!pt) return uissquare(u);
     308     2557659 :     if (!uissquareall(u, &a)) return 0;
     309     1181461 :     *pt = utoipos(a); return 1;
     310             :   }
     311     3361786 :   if (!squaremod(umodiu(x, 64*63*65*11))) return 0;
     312     2211486 :   av = avma; y = sqrtremi(x, &r);
     313     2211486 :   if (r != gen_0) return gc_long(av,0);
     314       40827 :   if (pt) { *pt = y; set_avma((pari_sp)y); } else set_avma(av);
     315       40827 :   return 1;
     316             : }
     317             : 
     318             : /* a t_INT, p prime */
     319             : long
     320      108937 : Zp_issquare(GEN a, GEN p)
     321             : {
     322             :   long v;
     323             :   GEN ap;
     324             : 
     325      108937 :   if (!signe(a) || equali1(a)) return 1;
     326      108937 :   v = Z_pvalrem(a, p, &ap);
     327      108937 :   if (v&1) return 0;
     328       61253 :   return absequaliu(p, 2)? Mod8(ap) == 1
     329       61253 :                          : kronecker(ap,p) == 1;
     330             : }
     331             : 
     332             : static long
     333        3710 : polissquareall(GEN x, GEN *pt)
     334             : {
     335             :   pari_sp av;
     336             :   long v;
     337             :   GEN y, a, b, p;
     338             : 
     339        3710 :   if (!signe(x))
     340             :   {
     341           7 :     if (pt) *pt = gcopy(x);
     342           7 :     return 1;
     343             :   }
     344        3703 :   if (odd(degpol(x))) return 0; /* odd degree */
     345        3703 :   av = avma;
     346        3703 :   v = RgX_valrem(x, &x);
     347        3703 :   if (v & 1) return gc_long(av,0);
     348        3696 :   a = gel(x,2); /* test constant coeff */
     349        3696 :   if (!pt)
     350          63 :   { if (!issquare(a)) return gc_long(av,0); }
     351             :   else
     352        3633 :   { if (!issquareall(a,&b)) return gc_long(av,0); }
     353        3696 :   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        3633 :   p = characteristic(x);
     358        3633 :   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        3598 :     long m = 1;
     382        3598 :     x = RgX_Rg_div(x,a);
     383             :     /* a(x^m) = B^2 => B = b(x^m) provided a(0) != 0 */
     384        3598 :     if (!signe(p)) x = RgX_deflate_max(x,&m);
     385        3598 :     y = ser2rfrac_i(gsqrt(RgX_to_ser(x,lg(x)-1),0));
     386        3605 :     if (!RgX_equal(RgX_sqr(y), x)) return gc_long(av,0);
     387        3591 :     if (!pt) return gc_long(av,1);
     388        3584 :     if (!gequal1(a)) y = gmul(b, y);
     389        3584 :     if (m != 1) y = RgX_inflate(y,m);
     390             :   }
     391        3626 : END:
     392        3626 :   if (v) y = RgX_shift_shallow(y, v>>1);
     393        3626 :   *pt = gerepilecopy(av, y); return 1;
     394             : }
     395             : 
     396             : /* b unit mod p */
     397             : static int
     398         637 : Up_ispower(GEN b, GEN K, GEN p, long d, GEN *pt)
     399             : {
     400         637 :   if (d == 1)
     401             :   { /* mod p: faster */
     402         553 :     if (!Fp_ispower(b, K, p)) return 0;
     403         553 :     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         616 :   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         777 : handle_pe(GEN *pa, GEN q, GEN L, GEN K, GEN p, long e)
     417             : {
     418             :   GEN t, A;
     419         777 :   long v = Z_pvalrem(*pa, p, &A), d = e - v;
     420         777 :   if (d <= 0) t = gen_0;
     421             :   else
     422             :   {
     423             :     ulong r;
     424         721 :     v = uabsdivui_rem(v, K, &r);
     425         721 :     if (r || !Up_ispower(A, K, p, d, L? &t: NULL)) return 0;
     426         616 :     if (L && v) t = mulii(t, powiu(p, v));
     427             :   }
     428         672 :   if (q) *pa = modii(*pa, q);
     429         672 :   if (L) vectrunc_append(L, mkintmod(t, powiu(p, e)));
     430         672 :   return 1;
     431             : }
     432             : long
     433         749 : 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         749 :   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         658 :   av = avma;
     451             : 
     452         658 :   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         658 :   if (!mod2(K)
     466         539 :       && kronecker(a, shifti(q,-vali(q))) == -1) return gc_long(av,0);
     467         651 :   L = pt? vectrunc_init(expi(q)+1): NULL;
     468         651 :   u_forprime_init(&S, 2, tridiv_bound(q));
     469      883911 :   while ((pp = u_forprime_next(&S)))
     470             :   {
     471             :     int stop;
     472      883757 :     e = Z_lvalrem_stop(&q, pp, &stop);
     473      883757 :     if (!e) continue;
     474         532 :     if (!handle_pe(&a, q, L, K, utoipos(pp), e)) return gc_long(av,0);
     475         511 :     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         154 :   l = lg(primetab);
     482         154 :   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         154 :   N = gcdii(a,q);
     491         154 :   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          84 :   if (!is_pm1(q))
     512             :   {
     513          84 :     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             :     {
     519          56 :       GEN part = ifac_start(q, 0);
     520             :       for(;;)
     521          84 :       {
     522             :         long e;
     523             :         GEN p;
     524         140 :         if (!ifac_next(&part, &p, &e)) break;
     525          98 :         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      251743 : issquareall(GEN x, GEN *pt)
     577             : {
     578      251743 :   long tx = typ(x);
     579             :   GEN F;
     580             :   pari_sp av;
     581             : 
     582      251743 :   if (!pt) return issquare(x);
     583       37568 :   switch(tx)
     584             :   {
     585       17058 :     case t_INT: return Z_issquareall(x, pt);
     586        1932 :     case t_FRAC: av = avma;
     587        1932 :       F = cgetg(3, t_FRAC);
     588        1932 :       if (   !Z_issquareall(gel(x,1), &gel(F,1))
     589        1932 :           || !Z_issquareall(gel(x,2), &gel(F,2))) return gc_long(av,0);
     590        1869 :       *pt = F; return 1;
     591             : 
     592          21 :     case t_POLMOD:
     593          21 :       return polmodispower(x, gen_2, pt);
     594        3640 :     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      229260 : issquare(GEN x)
     613             : {
     614             :   GEN a, p;
     615             :   long v;
     616             : 
     617      229260 :   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          77 :     case t_INTMOD:
     626          77 :       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         840 : polispower(GEN x, GEN K, GEN *pt)
     713             : {
     714             :   pari_sp av;
     715         840 :   long v, d, k = itos(K);
     716             :   GEN y, a, b;
     717         840 :   GEN T = NULL, p = NULL;
     718             : 
     719         840 :   if (!signe(x))
     720             :   {
     721           7 :     if (pt) *pt = gcopy(x);
     722           7 :     return 1;
     723             :   }
     724         833 :   d = degpol(x);
     725         833 :   if (d % k) return 0; /* degree not multiple of k */
     726         826 :   av = avma;
     727         826 :   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         490 :   v = RgX_valrem(x, &x);
     740         490 :   if (v % k) return 0;
     741         483 :   v /= k;
     742         483 :   a = gel(x,2); b = NULL;
     743         483 :   if (!ispower(a, K, &b)) return gc_long(av,0);
     744         469 :   if (d)
     745             :   {
     746         441 :     GEN p = characteristic(x);
     747         441 :     a = leading_coeff(x);
     748         441 :     if (!ispower(a, K, &b)) return gc_long(av,0);
     749         441 :     x = RgX_normalize(x);
     750         441 :     if (signe(p) && cmpii(p,K) <= 0)
     751           0 :       pari_err_IMPL("ispower(general t_POL) in small characteristic");
     752         441 :     y = gtrunc(gsqrtn(RgX_to_ser(x,lg(x)), K, NULL, 0));
     753         441 :     if (!RgX_equal(powgi(y, K), x)) return gc_long(av,0);
     754             :   }
     755             :   else
     756          28 :     y = pol_1(varn(x));
     757         469 :   if (pt)
     758             :   {
     759         434 :     if (!gequal1(a))
     760             :     {
     761          35 :       if (!b) b = gsqrtn(a, K, NULL, DEFAULTPREC);
     762          35 :       y = gmul(b,y);
     763             :     }
     764         434 :     if (v) y = RgX_shift_shallow(y, v);
     765         434 :     *pt = gerepilecopy(av, y);
     766             :   }
     767          35 :   else set_avma(av);
     768         469 :   return 1;
     769             : }
     770             : 
     771             : long
     772     1552287 : Z_ispowerall(GEN x, ulong k, GEN *pt)
     773             : {
     774     1552287 :   long s = signe(x);
     775             :   ulong mask;
     776     1552287 :   if (!s) { if (pt) *pt = gen_0; return 1; }
     777     1552287 :   if (s > 0) {
     778     1551958 :     if (k == 2) return Z_issquareall(x, pt);
     779     1423189 :     if (k == 3) { mask = 1; return !!is_357_power(x, pt, &mask); }
     780      227445 :     if (k == 5) { mask = 2; return !!is_357_power(x, pt, &mask); }
     781       59886 :     if (k == 7) { mask = 4; return !!is_357_power(x, pt, &mask); }
     782       53768 :     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         553 : Fp_ispower(GEN x, GEN K, GEN p)
     796             : {
     797         553 :   pari_sp av = avma;
     798             :   GEN p_1;
     799         553 :   x = modii(x, p);
     800         553 :   if (!signe(x) || equali1(x)) return gc_bool(av,1);
     801             :   /* implies p > 2 */
     802         112 :   p_1 = subiu(p,1);
     803         112 :   K = gcdii(K, p_1);
     804         112 :   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     8548893 : ispower(GEN x, GEN K, GEN *pt)
     853             : {
     854             :   GEN z;
     855             : 
     856     8548893 :   if (!K) return gisanypower(x, pt);
     857     1548711 :   if (typ(K) != t_INT) pari_err_TYPE("ispower",K);
     858     1548711 :   if (signe(K) <= 0) pari_err_DOMAIN("ispower","exponent","<=",gen_0,K);
     859     1548711 :   if (equali1(K)) { if (pt) *pt = gcopy(x); return 1; }
     860     1548662 :   switch(typ(x)) {
     861     1476035 :     case t_INT:
     862     1476035 :       if (lgefint(K) != 3) return 0;
     863     1475979 :       return Z_ispowerall(x, itou(K), pt);
     864       69981 :     case t_FRAC:
     865             :     {
     866       69981 :       GEN a = gel(x,1), b = gel(x,2);
     867             :       ulong k;
     868       69981 :       if (lgefint(K) != 3) return 0;
     869       69974 :       k = itou(K);
     870       69974 :       if (pt) {
     871       69911 :         z = cgetg(3, t_FRAC);
     872       69911 :         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       68427 :         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         840 :     case t_POL:
     889         840 :       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      505729 : split_exponent(ulong e, GEN *x)
     957             : {
     958             :   GEN fa, P, E;
     959      505729 :   long i, j, l, k = 1;
     960      505729 :   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             : static long
     978      865662 : Z_isanypower_nosmalldiv(GEN *px)
     979             : { /* any prime divisor of x is > 102 */
     980      865662 :   const double LOG2_103 = 6.6865; /* lower bound for log_2(103) */
     981      865662 :   const double LOG103 = 4.6347; /* lower bound for log(103) */
     982             :   forprime_t T;
     983      865662 :   ulong mask = 7, e2;
     984             :   long k, ex;
     985      865662 :   GEN y, x = *px;
     986             : 
     987      865662 :   k = 1;
     988      867664 :   while (Z_issquareall(x, &y)) { k <<= 1; x = y; }
     989      865882 :   while ( (ex = is_357_power(x, &y, &mask)) ) { k *= ex; x = y; }
     990      865662 :   e2 = (ulong)((expi(x) + 1) / LOG2_103); /* >= log_103 (x) */
     991      865662 :   if (u_forprime_init(&T, 11, e2))
     992             :   {
     993       17087 :     GEN logx = NULL;
     994       17087 :     const ulong Q = 30011; /* prime */
     995             :     ulong p, xmodQ;
     996       17087 :     double dlogx = 0;
     997             :     /* cut off at x^(1/p) ~ 2^30 bits which seems to be about optimum;
     998             :      * for large p the modular checks are no longer competitively fast */
     999       17129 :     while ( (ex = is_pth_power(x, &y, &T, 30)) )
    1000             :     {
    1001          42 :       k *= ex; x = y;
    1002          42 :       e2 = (ulong)((expi(x) + 1) / LOG2_103);
    1003          42 :       u_forprime_restrict(&T, e2);
    1004             :     }
    1005       17087 :     if (DEBUGLEVEL>4)
    1006           0 :       err_printf("Z_isanypower: now k=%ld, x=%ld-bit\n", k, expi(x)+1);
    1007       17087 :     xmodQ = umodiu(x, Q);
    1008             :     /* test Q | x, just in case */
    1009       17087 :     if (!xmodQ) { *px = x; return k * split_exponent(Z_lval(x,Q), px); }
    1010             :     /* x^(1/p) < 2^31 */
    1011       17073 :     p = T.p;
    1012       17073 :     if (p <= e2)
    1013             :     {
    1014       17059 :       logx = logr_abs( itor(x, DEFAULTPREC) );
    1015       17059 :       dlogx = rtodbl(logx);
    1016       17059 :       e2 = (ulong)(dlogx / LOG103); /* >= log_103(x) */
    1017             :     }
    1018      138530 :     while (p && p <= e2)
    1019             :     { /* is x a p-th power ? By computing y = round(x^(1/p)).
    1020             :        * Check whether y^p = x, first mod Q, then exactly. */
    1021      121457 :       pari_sp av = avma;
    1022             :       long e;
    1023      121457 :       GEN logy = divru(logx, p), y = grndtoi(mpexp(logy), &e);
    1024      121457 :       ulong ymodQ = umodiu(y,Q);
    1025      121457 :       if (e >= -10 || Fl_powu(ymodQ, p % (Q-1), Q) != xmodQ
    1026      121457 :                    || !equalii(powiu(y, p), x)) set_avma(av);
    1027             :       else
    1028             :       {
    1029          21 :         k *= p; x = y; xmodQ = ymodQ; logx = logy; dlogx /= p;
    1030          21 :         e2 = (ulong)(dlogx / LOG103); /* >= log_103(x) */
    1031          21 :         u_forprime_restrict(&T, e2);
    1032          21 :         continue; /* if success, retry same p */
    1033             :       }
    1034      121436 :       p = u_forprime_next(&T);
    1035             :     }
    1036             :   }
    1037      865648 :   *px = x; return k;
    1038             : }
    1039             : 
    1040             : static ulong tinyprimes[] = {
    1041             :   2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71,
    1042             :   73, 79, 83, 89, 97, 101, 103, 107, 109, 113, 127, 131, 137, 139, 149, 151,
    1043             :   157, 163, 167, 173, 179, 181, 191, 193, 197, 199
    1044             : };
    1045             : 
    1046             : /* disregard the sign of x, caller will take care of x < 0 */
    1047             : static long
    1048     7001890 : Z_isanypower_aux(GEN x, GEN *pty)
    1049             : {
    1050             :   long ex, v, i, l, k;
    1051             :   GEN y, P, E;
    1052     7001890 :   ulong mask, e = 0;
    1053             : 
    1054     7001890 :   if (abscmpii(x, gen_2) < 0) return 0; /* -1,0,1 */
    1055             : 
    1056     7001876 :   if (signe(x) < 0) x = negi(x);
    1057     7001876 :   k = l = 1;
    1058     7001876 :   P = cgetg(26 + 1, t_VECSMALL);
    1059     7001876 :   E = cgetg(26 + 1, t_VECSMALL);
    1060             :   /* trial division */
    1061    61498458 :   for(i = 0; i < 26; i++)
    1062             :   {
    1063    60162963 :     ulong p = tinyprimes[i];
    1064             :     int stop;
    1065    60162963 :     v = Z_lvalrem_stop(&x, p, &stop);
    1066    60162963 :     if (v)
    1067             :     {
    1068     7922376 :       P[l] = p;
    1069     7922376 :       E[l] = v; l++;
    1070     8170456 :       e = ugcd(e, v); if (e == 1) goto END;
    1071             :     }
    1072    54744662 :     if (stop) {
    1073      248080 :       if (is_pm1(x)) k = e;
    1074      248080 :       goto END;
    1075             :     }
    1076             :   }
    1077             : 
    1078     1335495 :   if (e)
    1079             :   { /* Bingo. Result divides e */
    1080             :     long v3, v5, v7;
    1081      505715 :     ulong e2 = e;
    1082      505715 :     v = u_lvalrem(e2, 2, &e2);
    1083      505715 :     if (v)
    1084             :     {
    1085      375277 :       for (i = 0; i < v; i++)
    1086             :       {
    1087      374185 :         if (!Z_issquareall(x, &y)) break;
    1088        1302 :         k <<= 1; x = y;
    1089             :       }
    1090             :     }
    1091      505715 :     mask = 0;
    1092      505715 :     v3 = u_lvalrem(e2, 3, &e2); if (v3) mask = 1;
    1093      505715 :     v5 = u_lvalrem(e2, 5, &e2); if (v5) mask |= 2;
    1094      505715 :     v7 = u_lvalrem(e2, 7, &e2); if (v7) mask |= 4;
    1095     1011507 :     while ( (ex = is_357_power(x, &y, &mask)) ) {
    1096          77 :       x = y;
    1097          77 :       switch(ex)
    1098             :       {
    1099          28 :         case 3: k *= 3; if (--v3 == 0) mask &= ~1; break;
    1100          28 :         case 5: k *= 5; if (--v5 == 0) mask &= ~2; break;
    1101          21 :         case 7: k *= 7; if (--v7 == 0) mask &= ~4; break;
    1102             :       }
    1103      505792 :     }
    1104      505715 :     k *= split_exponent(e2, &x);
    1105             :   }
    1106             :   else
    1107      829780 :     k = Z_isanypower_nosmalldiv(&x);
    1108     7001876 : END:
    1109     7001876 :   if (pty && k != 1)
    1110             :   {
    1111        8743 :     if (e)
    1112             :     { /* add missing small factors */
    1113        6881 :       y = powuu(P[1], E[1] / k);
    1114       14035 :       for (i = 2; i < l; i++) y = mulii(y, powuu(P[i], E[i] / k));
    1115        6881 :       x = equali1(x)? y: mulii(x,y);
    1116             :     }
    1117        8743 :     *pty = x;
    1118             :   }
    1119     7001876 :   return k == 1? 0: k;
    1120             : }
    1121             : 
    1122             : long
    1123     7001890 : Z_isanypower(GEN x, GEN *pty)
    1124             : {
    1125     7001890 :   pari_sp av = avma;
    1126     7001890 :   long k = Z_isanypower_aux(x, pty);
    1127     7001890 :   if (!k) return gc_long(av,0);
    1128        8806 :   if (signe(x) < 0)
    1129             :   {
    1130          42 :     long v = vals(k);
    1131          42 :     if (v)
    1132             :     {
    1133          28 :       k >>= v;
    1134          28 :       if (k == 1) return gc_long(av,0);
    1135          21 :       if (!pty) return gc_long(av,k);
    1136          14 :       *pty = gerepileuptoint(av, powiu(*pty, 1<<v));
    1137          14 :       togglesign(*pty); return k;
    1138             :     }
    1139          14 :     if (pty) togglesign_safe(pty);
    1140             :   }
    1141        8778 :   if (!pty) return gc_long(av, k);
    1142        8722 :   *pty = gerepilecopy(av, *pty); return k;
    1143             : }
    1144             : 
    1145             : /* Faster than expi(n) == vali(n) or hamming(n) == 1 even for single-word
    1146             :  * values. If all you have is a word, you can just use n & !(n & (n-1)). */
    1147             : long
    1148      529599 : Z_ispow2(GEN n)
    1149             : {
    1150             :   GEN xp;
    1151             :   long i, l;
    1152             :   ulong u;
    1153      529599 :   if (signe(n) != 1) return 0;
    1154      473592 :   xp = int_LSW(n); u = *xp; l = lgefint(n);
    1155      534473 :   for (i = 3; i < l; ++i)
    1156             :   {
    1157      253551 :     if (u) return 0;
    1158       60881 :     xp = int_nextW(xp); u = *xp;
    1159             :   }
    1160      280922 :   return !(u & (u-1));
    1161             : }
    1162             : 
    1163             : static long
    1164      842120 : isprimepower_i(GEN n, GEN *pt, long flag)
    1165             : {
    1166      842120 :   pari_sp av = avma;
    1167             :   long i, v;
    1168             : 
    1169      842120 :   if (typ(n) != t_INT) pari_err_TYPE("isprimepower", n);
    1170      842120 :   if (signe(n) <= 0) return 0;
    1171             : 
    1172      842120 :   if (lgefint(n) == 3)
    1173             :   {
    1174             :     ulong p;
    1175      541225 :     v = uisprimepower(n[2], &p);
    1176      541225 :     if (v)
    1177             :     {
    1178       54999 :       if (pt) *pt = utoipos(p);
    1179       54999 :       return v;
    1180             :     }
    1181      486226 :     return 0;
    1182             :   }
    1183     1663672 :   for (i = 0; i < 26; i++)
    1184             :   {
    1185     1627790 :     ulong p = tinyprimes[i];
    1186     1627790 :     v = Z_lvalrem(n, p, &n);
    1187     1627790 :     if (v)
    1188             :     {
    1189      265013 :       set_avma(av);
    1190      265013 :       if (!is_pm1(n)) return 0;
    1191         617 :       if (pt) *pt = utoipos(p);
    1192         617 :       return v;
    1193             :     }
    1194             :   }
    1195             :   /* p | n => p >= 103 */
    1196       35882 :   v = Z_isanypower_nosmalldiv(&n); /* expensive */
    1197       35882 :   if (!(flag? isprime(n): BPSW_psp(n))) return gc_long(av,0);
    1198        5570 :   if (pt) *pt = gerepilecopy(av, n); else set_avma(av);
    1199        5570 :   return v;
    1200             : }
    1201             : long
    1202      840147 : isprimepower(GEN n, GEN *pt) { return isprimepower_i(n,pt,1); }
    1203             : long
    1204        1973 : ispseudoprimepower(GEN n, GEN *pt) { return isprimepower_i(n,pt,0); }
    1205             : 
    1206             : long
    1207      642137 : uisprimepower(ulong n, ulong *pp)
    1208             : { /* We must have CUTOFF^11 >= ULONG_MAX and CUTOFF^3 < ULONG_MAX.
    1209             :    * Tests suggest that 200-300 is the best range for 64-bit platforms. */
    1210      642137 :   const ulong CUTOFF = 200UL;
    1211      642137 :   const long TINYCUTOFF = 46;  /* tinyprimes[45] = 199 */
    1212      642137 :   const ulong CUTOFF3 = CUTOFF*CUTOFF*CUTOFF;
    1213             : #ifdef LONG_IS_64BIT
    1214             :   /* primes preceeding the appropriate root of ULONG_MAX. */
    1215      567558 :   const ulong ROOT9 = 137;
    1216      567558 :   const ulong ROOT8 = 251;
    1217      567558 :   const ulong ROOT7 = 563;
    1218      567558 :   const ulong ROOT5 = 7129;
    1219      567558 :   const ulong ROOT4 = 65521;
    1220             : #else
    1221       74579 :   const ulong ROOT9 = 11;
    1222       74579 :   const ulong ROOT8 = 13;
    1223       74579 :   const ulong ROOT7 = 23;
    1224       74579 :   const ulong ROOT5 = 83;
    1225       74579 :   const ulong ROOT4 = 251;
    1226             : #endif
    1227             :   ulong mask;
    1228             :   long v, i;
    1229             :   int e;
    1230      642137 :   if (n < 2) return 0;
    1231      642102 :   if (!odd(n)) {
    1232      341376 :     if (n & (n-1)) return 0;
    1233       65111 :     *pp = 2; return vals(n);
    1234             :   }
    1235      300726 :   if (n < 8) { *pp = n; return 1; } /* 3,5,7 */
    1236     3655068 :   for (i = 1/*skip p=2*/; i < TINYCUTOFF; i++)
    1237             :   {
    1238     3595993 :     ulong p = tinyprimes[i];
    1239     3595993 :     if (n % p == 0)
    1240             :     {
    1241      212426 :       v = u_lvalrem(n, p, &n);
    1242      212426 :       if (n == 1) { *pp = p; return v; }
    1243      209998 :       return 0;
    1244             :     }
    1245             :   }
    1246             :   /* p | n => p >= CUTOFF */
    1247             : 
    1248       59075 :   if (n < CUTOFF3)
    1249             :   {
    1250       46354 :     if (n < CUTOFF*CUTOFF || uisprime_101(n)) { *pp = n; return 1; }
    1251           0 :     if (uissquareall(n, &n)) { *pp = n; return 2; }
    1252           0 :     return 0;
    1253             :   }
    1254             : 
    1255             :   /* Check for squares, fourth powers, and eighth powers as appropriate. */
    1256       12721 :   v = 1;
    1257       12721 :   if (uissquareall(n, &n)) {
    1258           0 :     v <<= 1;
    1259           0 :     if (CUTOFF <= ROOT4 && uissquareall(n, &n)) {
    1260           0 :       v <<= 1;
    1261           0 :       if (CUTOFF <= ROOT8 && uissquareall(n, &n)) v <<= 1;
    1262             :     }
    1263             :   }
    1264             : 
    1265       12721 :   if (CUTOFF > ROOT5) mask = 1;
    1266             :   else
    1267             :   {
    1268       12720 :     const ulong CUTOFF5 = CUTOFF3*CUTOFF*CUTOFF;
    1269       12720 :     if (n < CUTOFF5) mask = 1; else mask = 3;
    1270       12720 :     if (CUTOFF <= ROOT7)
    1271             :     {
    1272       12720 :       const ulong CUTOFF7 = CUTOFF5*CUTOFF*CUTOFF;
    1273       12720 :       if (n >= CUTOFF7) mask = 7;
    1274             :     }
    1275             :   }
    1276             : 
    1277       12721 :   if (CUTOFF <= ROOT9 && (e = uis_357_power(n, &n, &mask))) { v *= e; mask=1; }
    1278       12721 :   if ((e = uis_357_power(n, &n, &mask))) v *= e;
    1279             : 
    1280       12721 :   if (uisprime_101(n)) { *pp = n; return v; }
    1281        6984 :   return 0;
    1282             : }
    1283             : 

Generated by: LCOV version 1.14