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 29395-ef22f77854) Lines: 685 735 93.2 %
Date: 2024-06-15 09:03:40 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      829329 : ulogintall(ulong B, ulong y, ulong *ptq)
      29             : {
      30             :   ulong r, r2;
      31             :   long e;
      32             : 
      33      829329 :   if (y == 2)
      34             :   {
      35       22892 :     long eB = expu(B); /* 2^eB <= B < 2^(eB + 1) */
      36       22892 :     if (ptq) *ptq = 1UL << eB;
      37       22892 :     return eB;
      38             :   }
      39      806437 :   r = y, r2 = 1UL;
      40     2863996 :   for (e=1;; e++)
      41             :   { /* here, r = y^e, r2 = y^(e-1) */
      42     2863996 :     if (r >= B)
      43             :     {
      44      804619 :       if (r != B) { e--; r = r2; }
      45      804619 :       if (ptq) *ptq = r;
      46      804619 :       return e;
      47             :     }
      48     2059377 :     r2 = r;
      49     2059377 :     r = umuluu_or_0(y, r);
      50     2059386 :     if (!r)
      51             :     {
      52        1827 :       if (ptq) *ptq = r2;
      53        1827 :       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      928965 : logintall(GEN B, GEN y, GEN *ptq)
      62             : {
      63             :   pari_sp av;
      64      928965 :   long ey, e, emax, i, eB = expi(B); /* 2^eB <= B < 2^(eB + 1) */
      65             :   GEN q, pow2;
      66             : 
      67      928963 :   if (lgefint(B) == 3)
      68             :   {
      69             :     ulong q;
      70      828180 :     if (lgefint(y) > 3)
      71             :     {
      72           0 :       if (ptq) *ptq = gen_1;
      73           0 :       return 0;
      74             :     }
      75      828180 :     if (!ptq) return ulogintall(B[2], y[2], NULL);
      76      126475 :     e = ulogintall(B[2], y[2], &q);
      77      126477 :     *ptq = utoi(q); return e;
      78             :   }
      79      100783 :   if (equaliu(y,2))
      80             :   {
      81         855 :     if (ptq) *ptq = int2n(eB);
      82         855 :     return eB;
      83             :   }
      84       99942 :   av = avma;
      85       99942 :   ey = expi(y);
      86             :   /* eB/(ey+1) - 1 < e <= eB/ey */
      87       99941 :   emax = eB/ey;
      88       99941 :   if (emax <= 13) /* e small, be naive */
      89             :   {
      90       10842 :     GEN r = y, r2 = gen_1;
      91       10842 :     for (e=1;; e++)
      92      104891 :     { /* here, r = y^e, r2 = y^(e-1) */
      93      115733 :       long fl = cmpii(r, B);
      94      115733 :       if (fl >= 0)
      95             :       {
      96       10842 :         if (fl) { e--; cgiv(r); r = r2; }
      97       10842 :         if (ptq) *ptq = gerepileuptoint(av, r); else set_avma(av);
      98       10842 :         return e;
      99             :       }
     100      104891 :       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       89099 :   pow2 = new_chunk((long)log2(eB)+2);
     108       89099 :   gel(pow2,0) = y;
     109       89099 :   for (i=0, q=y;; )
     110      416992 :   {
     111      506091 :     GEN r = gel(pow2,i); /* r = y^2^i */
     112      506091 :     long fl = cmpii(r,B);
     113      506065 :     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      506065 :     if (fl > 0) { i--; break; }
     120      474587 :     q = r;
     121      474587 :     if (1L<<(i+1) > emax) break;
     122      416962 :     gel(pow2,++i) = sqri(q);
     123             :   }
     124             : 
     125       89103 :   for (e = 1L<<i;;)
     126      385514 :   { /* y^e = q < B < r = q * y^(2^i) */
     127      474617 :     pari_sp av2 = avma;
     128             :     long fl;
     129             :     GEN r;
     130      474617 :     if (--i < 0) break;
     131      385523 :     r = mulii(q, gel(pow2,i));
     132      385530 :     fl = cmpii(r, B);
     133      385553 :     if (fl > 0) set_avma(av2);
     134             :     else
     135             :     {
     136      174048 :       e += (1L<<i);
     137      174048 :       q = r;
     138      174048 :       if (!fl) break; /* B = r */
     139             :     }
     140             :   }
     141       89101 :   if (ptq) *ptq = gerepileuptoint(av, q); else set_avma(av);
     142       89101 :   return e;
     143             : }
     144             : 
     145             : long
     146       26601 : logint0(GEN B, GEN y, GEN *ptq)
     147             : {
     148       26601 :   const char *f = "logint";
     149       26601 :   if (typ(y) != t_INT) pari_err_TYPE(f,y);
     150       26601 :   if (cmpis(y, 2) < 0) pari_err_DOMAIN(f, "b" ,"<=", gen_1, y);
     151       26601 :   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       26524 :   if (signe(B) <= 0) pari_err_DOMAIN(f, "x" ,"<=", gen_0, B);
     189       26524 :   return logintall(B,y,ptq);
     190             : }
     191             : 
     192             : /*********************************************************************/
     193             : /**                     INTEGRAL SQUARE ROOT                        **/
     194             : /*********************************************************************/
     195             : GEN
     196       88023 : sqrtint(GEN a)
     197             : {
     198       88023 :   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       87967 :   switch (signe(a))
     225             :   {
     226       87946 :     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    14101841 : squaremod(ulong A)
     256             : {
     257    14101841 :   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    14101841 :   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    14101841 :   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    14101841 :   const int squaremod11[]={1,1,0,1,1,1,0,0,0,1, 0};
     267    14101841 :   return (squaremod64[A & 0x3fUL]
     268     6707055 :     && squaremod63[A % 63UL]
     269     4689776 :     && squaremod65[A % 65UL]
     270    20808896 :     && squaremod11[A % 11UL]);
     271             : }
     272             : 
     273             : /* emulate Z_issquareall on single-word integers */
     274             : long
     275    10465530 : uissquareall(ulong A, ulong *sqrtA)
     276             : {
     277    10465530 :   if (!A) { *sqrtA = 0; return 1; }
     278    10465530 :   if (squaremod(A))
     279             :   {
     280     1932216 :     ulong a = usqrt(A);
     281     1932204 :     if (a * a == A) { *sqrtA = a; return 1; }
     282             :   }
     283     8562367 :   return 0;
     284             : }
     285             : long
     286      274315 : uissquare(ulong A)
     287             : {
     288      274315 :   if (!A) return 1;
     289      274315 :   if (squaremod(A)) { ulong a = usqrt(A); if (a * a == A) return 1; }
     290      264571 :   return 0;
     291             : }
     292             : 
     293             : long
     294     6678334 : Z_issquareall(GEN x, GEN *pt)
     295             : {
     296             :   pari_sp av;
     297             :   GEN y, r;
     298             : 
     299     6678334 :   switch(signe(x))
     300             :   {
     301      180654 :     case -1: return 0;
     302        1806 :     case 0: if (pt) *pt=gen_0; return 1;
     303             :   }
     304     6495874 :   if (lgefint(x) == 3)
     305             :   {
     306     3133885 :     ulong u = uel(x,2), a;
     307     3133885 :     if (!pt) return uissquare(u);
     308     2859570 :     if (!uissquareall(u, &a)) return 0;
     309     1476133 :     *pt = utoipos(a); return 1;
     310             :   }
     311     3361989 :   if (!squaremod(umodiu(x, 64*63*65*11))) return 0;
     312     2211501 :   av = avma; y = sqrtremi(x, &r);
     313     2211501 :   if (r != gen_0) return gc_long(av,0);
     314       40951 :   if (pt) { *pt = y; set_avma((pari_sp)y); } else set_avma(av);
     315       40951 :   return 1;
     316             : }
     317             : 
     318             : /* a t_INT, p prime */
     319             : long
     320      109728 : Zp_issquare(GEN a, GEN p)
     321             : {
     322             :   long v;
     323             :   GEN ap;
     324             : 
     325      109728 :   if (!signe(a) || equali1(a)) return 1;
     326      109728 :   v = Z_pvalrem(a, p, &ap);
     327      109728 :   if (v&1) return 0;
     328       62114 :   return absequaliu(p, 2)? Mod8(ap) == 1
     329       62114 :                          : 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     4104231 :   while ((pp = u_forprime_next(&S)))
     470             :   {
     471             :     int stop;
     472     4104077 :     e = Z_lvalrem_stop(&q, pp, &stop);
     473     4104077 :     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      251771 : issquareall(GEN x, GEN *pt)
     577             : {
     578      251771 :   long tx = typ(x);
     579             :   GEN F;
     580             :   pari_sp av;
     581             : 
     582      251771 :   if (!pt) return issquare(x);
     583       37596 :   switch(tx)
     584             :   {
     585       17079 :     case t_INT: return Z_issquareall(x, pt);
     586        1939 :     case t_FRAC: av = avma;
     587        1939 :       F = cgetg(3, t_FRAC);
     588        1939 :       if (   !Z_issquareall(gel(x,1), &gel(F,1))
     589        1939 :           || !Z_issquareall(gel(x,2), &gel(F,2))) return gc_long(av,0);
     590        1876 :       *pt = F; return 1;
     591             : 
     592          21 :     case t_POLMOD:
     593          21 :       return polmodispower(x, gen_2, pt);
     594        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     1558551 : Z_ispowerall(GEN x, ulong k, GEN *pt)
     773             : {
     774     1558551 :   long s = signe(x);
     775             :   ulong mask;
     776     1558551 :   if (!s) { if (pt) *pt = gen_0; return 1; }
     777     1558551 :   if (s > 0) {
     778     1558222 :     if (k == 2) return Z_issquareall(x, pt);
     779     1425197 :     if (k == 3) { mask = 1; return !!is_357_power(x, pt, &mask); }
     780      227493 :     if (k == 5) { mask = 2; return !!is_357_power(x, pt, &mask); }
     781       59920 :     if (k == 7) { mask = 4; return !!is_357_power(x, pt, &mask); }
     782       53802 :     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     8555178 : ispower(GEN x, GEN K, GEN *pt)
     853             : {
     854             :   GEN z;
     855             : 
     856     8555178 :   if (!K) return gisanypower(x, pt);
     857     1554996 :   if (typ(K) != t_INT) pari_err_TYPE("ispower",K);
     858     1554996 :   if (signe(K) <= 0) pari_err_DOMAIN("ispower","exponent","<=",gen_0,K);
     859     1554996 :   if (equali1(K)) { if (pt) *pt = gcopy(x); return 1; }
     860     1554947 :   switch(typ(x)) {
     861     1482320 :     case t_INT:
     862     1482320 :       if (lgefint(K) != 3) return 0;
     863     1482264 :       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      507453 : split_exponent(ulong e, GEN *x)
     957             : {
     958             :   GEN fa, P, E;
     959      507453 :   long i, j, l, k = 1;
     960      507453 :   if (e == 1) return 1;
     961          14 :   fa = factoru(e);
     962          14 :   P = gel(fa,1);
     963          14 :   E = gel(fa,2); l = lg(P);
     964          28 :   for (i = 1; i < l; i++)
     965             :   {
     966          14 :     ulong p = P[i];
     967          28 :     for (j = 0; j < E[i]; j++)
     968             :     {
     969             :       GEN y;
     970          14 :       if (!is_kth_power(*x, p, &y)) break;
     971          14 :       k *= p; *x = y;
     972             :     }
     973             :   }
     974          14 :   return k;
     975             : }
     976             : 
     977             : /* any prime divisor of x is > 102 */
     978             : static long
     979      869930 : Z_isanypower_101(GEN *px)
     980             : {
     981      869930 :   const double LOG2_103 = 6.6865; /* lower bound for log_2(103) */
     982      869930 :   const double LOG103 = 4.6347; /* lower bound for log(103) */
     983             :   forprime_t T;
     984      869930 :   ulong mask = 7, e2;
     985             :   long k, ex;
     986      869930 :   GEN y, x = *px;
     987             : 
     988      869930 :   k = 1;
     989      873135 :   while (Z_issquareall(x, &y)) { k <<= 1; x = y; }
     990      870358 :   while ( (ex = is_357_power(x, &y, &mask)) ) { k *= ex; x = y; }
     991      869930 :   e2 = (ulong)((expi(x) + 1) / LOG2_103); /* >= log_103 (x) */
     992      869930 :   if (u_forprime_init(&T, 11, e2))
     993             :   {
     994       17381 :     GEN logx = NULL;
     995       17381 :     const ulong Q = 30011; /* prime */
     996             :     ulong p, xmodQ;
     997       17381 :     double dlogx = 0;
     998             :     /* cut off at x^(1/p) ~ 2^30 bits which seems to be about optimum;
     999             :      * for large p the modular checks are no longer competitively fast */
    1000       17423 :     while ( (ex = is_pth_power(x, &y, &T, 30)) )
    1001             :     {
    1002          42 :       k *= ex; x = y;
    1003          42 :       e2 = (ulong)((expi(x) + 1) / LOG2_103);
    1004          42 :       u_forprime_restrict(&T, e2);
    1005             :     }
    1006       17381 :     if (DEBUGLEVEL>4)
    1007           0 :       err_printf("Z_isanypower: now k=%ld, x=%ld-bit\n", k, expi(x)+1);
    1008       17381 :     xmodQ = umodiu(x, Q);
    1009             :     /* test Q | x, just in case */
    1010       17381 :     if (!xmodQ) { *px = x; return k * split_exponent(Z_lval(x,Q), px); }
    1011             :     /* x^(1/p) < 2^31 */
    1012       17367 :     p = T.p;
    1013       17367 :     if (p <= e2)
    1014             :     {
    1015       17353 :       logx = logr_abs( itor(x, DEFAULTPREC) );
    1016       17353 :       dlogx = rtodbl(logx);
    1017       17353 :       e2 = (ulong)(dlogx / LOG103); /* >= log_103(x) */
    1018             :     }
    1019      143171 :     while (p && p <= e2)
    1020             :     { /* is x a p-th power ? By computing y = round(x^(1/p)).
    1021             :        * Check whether y^p = x, first mod Q, then exactly. */
    1022      125804 :       pari_sp av = avma;
    1023             :       long e;
    1024      125804 :       GEN logy = divru(logx, p), y = grndtoi(mpexp(logy), &e);
    1025      125804 :       ulong ymodQ = umodiu(y,Q);
    1026      125804 :       if (e >= -10 || Fl_powu(ymodQ, p % (Q-1), Q) != xmodQ
    1027      125804 :                    || !equalii(powiu(y, p), x)) set_avma(av);
    1028             :       else
    1029             :       {
    1030          21 :         k *= p; x = y; xmodQ = ymodQ; logx = logy; dlogx /= p;
    1031          21 :         e2 = (ulong)(dlogx / LOG103); /* >= log_103(x) */
    1032          21 :         u_forprime_restrict(&T, e2);
    1033          21 :         continue; /* if success, retry same p */
    1034             :       }
    1035      125783 :       p = u_forprime_next(&T);
    1036             :     }
    1037             :   }
    1038      869916 :   *px = x; return k;
    1039             : }
    1040             : 
    1041             : static ulong tinyprimes[] = {
    1042             :   2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71,
    1043             :   73, 79, 83, 89, 97, 101, 103, 107, 109, 113, 127, 131, 137, 139, 149, 151,
    1044             :   157, 163, 167, 173, 179, 181, 191, 193, 197, 199
    1045             : };
    1046             : 
    1047             : /* disregard the sign of x, caller will take care of x < 0 */
    1048             : static long
    1049     7195282 : Z_isanypower_aux(GEN x, GEN *pty)
    1050             : {
    1051             :   long ex, v, i, l, k;
    1052             :   GEN y, P, E;
    1053     7195282 :   ulong mask, e = 0;
    1054             : 
    1055     7195282 :   if (abscmpii(x, gen_2) < 0) return 0; /* -1,0,1 */
    1056             : 
    1057     7194744 :   x = absi_shallow(x);
    1058     7194744 :   k = l = 1;
    1059     7194744 :   P = cgetg(26 + 1, t_VECSMALL);
    1060     7194744 :   E = cgetg(26 + 1, t_VECSMALL);
    1061             :   /* trial division */
    1062    61981272 :   for(i = 0; i < 26; i++)
    1063             :   {
    1064    60639785 :     ulong p = tinyprimes[i];
    1065             :     int stop;
    1066    60639785 :     v = Z_lvalrem_stop(&x, p, &stop);
    1067    60639785 :     if (v)
    1068             :     {
    1069     8043920 :       P[l] = p;
    1070     8043920 :       E[l] = v; l++;
    1071     8423896 :       e = ugcd(e, v); if (e == 1) goto END;
    1072             :     }
    1073    55166504 :     if (stop) {
    1074      379976 :       if (is_pm1(x)) k = e;
    1075      379976 :       goto END;
    1076             :     }
    1077             :   }
    1078             : 
    1079     1341487 :   if (e)
    1080             :   { /* Bingo. Result divides e */
    1081             :     long v3, v5, v7;
    1082      507439 :     ulong e2 = e;
    1083      507439 :     v = u_lvalrem(e2, 2, &e2);
    1084      507439 :     if (v)
    1085             :     {
    1086      377776 :       for (i = 0; i < v; i++)
    1087             :       {
    1088      375661 :         if (!Z_issquareall(x, &y)) break;
    1089        2331 :         k <<= 1; x = y;
    1090             :       }
    1091             :     }
    1092      507439 :     mask = 0;
    1093      507439 :     v3 = u_lvalrem(e2, 3, &e2); if (v3) mask = 1;
    1094      507439 :     v5 = u_lvalrem(e2, 5, &e2); if (v5) mask |= 2;
    1095      507439 :     v7 = u_lvalrem(e2, 7, &e2); if (v7) mask |= 4;
    1096      507522 :     while ( (ex = is_357_power(x, &y, &mask)) ) {
    1097          83 :       x = y;
    1098          83 :       switch(ex)
    1099             :       {
    1100          34 :         case 3: k *= 3; if (--v3 == 0) mask &= ~1; break;
    1101          28 :         case 5: k *= 5; if (--v5 == 0) mask &= ~2; break;
    1102          21 :         case 7: k *= 7; if (--v7 == 0) mask &= ~4; break;
    1103             :       }
    1104             :     }
    1105      507439 :     k *= split_exponent(e2, &x);
    1106             :   }
    1107             :   else
    1108      834048 :     k = Z_isanypower_101(&x);
    1109     7194744 : END:
    1110     7194744 :   if (pty && k != 1)
    1111             :   {
    1112       69965 :     if (e)
    1113             :     { /* add missing small factors */
    1114       66879 :       y = powuu(P[1], E[1] / k);
    1115       76877 :       for (i = 2; i < l; i++) y = mulii(y, powuu(P[i], E[i] / k));
    1116       66879 :       x = equali1(x)? y: mulii(x,y);
    1117             :     }
    1118       69965 :     *pty = x;
    1119             :   }
    1120     7194744 :   return k == 1? 0: k;
    1121             : }
    1122             : 
    1123             : long
    1124     7195282 : Z_isanypower(GEN x, GEN *pty)
    1125             : {
    1126     7195282 :   pari_sp av = avma;
    1127     7195282 :   long k = Z_isanypower_aux(x, pty);
    1128     7195282 :   if (!k) return gc_long(av,0);
    1129       70028 :   if (signe(x) < 0)
    1130             :   {
    1131          42 :     long v = vals(k);
    1132          42 :     if (v)
    1133             :     {
    1134          28 :       k >>= v;
    1135          28 :       if (k == 1) return gc_long(av,0);
    1136          21 :       if (!pty) return gc_long(av,k);
    1137          14 :       *pty = gerepileuptoint(av, powiu(*pty, 1<<v));
    1138          14 :       togglesign(*pty); return k;
    1139             :     }
    1140          14 :     if (pty) togglesign_safe(pty);
    1141             :   }
    1142       70000 :   if (!pty) return gc_long(av, k);
    1143       69944 :   *pty = gerepilecopy(av, *pty); return k;
    1144             : }
    1145             : 
    1146             : /* Faster than expi(n) == vali(n) or hamming(n) == 1 even for single-word
    1147             :  * values. If all you have is a word, you can just use n & !(n & (n-1)). */
    1148             : long
    1149      529706 : Z_ispow2(GEN n)
    1150             : {
    1151             :   GEN xp;
    1152             :   long i, l;
    1153             :   ulong u;
    1154      529706 :   if (signe(n) != 1) return 0;
    1155      473699 :   xp = int_LSW(n); u = *xp; l = lgefint(n);
    1156      534580 :   for (i = 3; i < l; ++i)
    1157             :   {
    1158      253658 :     if (u) return 0;
    1159       60881 :     xp = int_nextW(xp); u = *xp;
    1160             :   }
    1161      280922 :   return !(u & (u-1));
    1162             : }
    1163             : 
    1164             : static long
    1165      842128 : isprimepower_i(GEN n, GEN *pt, long flag)
    1166             : {
    1167      842128 :   pari_sp av = avma;
    1168             :   long i, v;
    1169             : 
    1170      842128 :   if (typ(n) != t_INT) pari_err_TYPE("isprimepower", n);
    1171      842128 :   if (signe(n) <= 0) return 0;
    1172             : 
    1173      842128 :   if (lgefint(n) == 3)
    1174             :   {
    1175             :     ulong p;
    1176      541225 :     v = uisprimepower(n[2], &p);
    1177      541225 :     if (v)
    1178             :     {
    1179       54999 :       if (pt) *pt = utoipos(p);
    1180       54999 :       return v;
    1181             :     }
    1182      486226 :     return 0;
    1183             :   }
    1184     1663680 :   for (i = 0; i < 26; i++)
    1185             :   {
    1186     1627798 :     ulong p = tinyprimes[i];
    1187     1627798 :     v = Z_lvalrem(n, p, &n);
    1188     1627798 :     if (v)
    1189             :     {
    1190      265021 :       set_avma(av);
    1191      265021 :       if (!is_pm1(n)) return 0;
    1192         625 :       if (pt) *pt = utoipos(p);
    1193         625 :       return v;
    1194             :     }
    1195             :   }
    1196             :   /* p | n => p >= 103 */
    1197       35882 :   v = Z_isanypower_101(&n); /* expensive */
    1198       35882 :   if (!(flag? isprime(n): BPSW_psp(n))) return gc_long(av,0);
    1199        5570 :   if (pt) *pt = gerepilecopy(av, n); else set_avma(av);
    1200        5570 :   return v;
    1201             : }
    1202             : long
    1203      840147 : isprimepower(GEN n, GEN *pt) { return isprimepower_i(n,pt,1); }
    1204             : long
    1205        1981 : ispseudoprimepower(GEN n, GEN *pt) { return isprimepower_i(n,pt,0); }
    1206             : 
    1207             : long
    1208      642176 : uisprimepower(ulong n, ulong *pp)
    1209             : { /* We must have CUTOFF^11 >= ULONG_MAX and CUTOFF^3 < ULONG_MAX.
    1210             :    * Tests suggest that 200-300 is the best range for 64-bit platforms. */
    1211      642176 :   const ulong CUTOFF = 200UL;
    1212      642176 :   const long TINYCUTOFF = 46;  /* tinyprimes[45] = 199 */
    1213      642176 :   const ulong CUTOFF3 = CUTOFF*CUTOFF*CUTOFF;
    1214             : #ifdef LONG_IS_64BIT
    1215             :   /* primes preceeding the appropriate root of ULONG_MAX. */
    1216      567591 :   const ulong ROOT9 = 137;
    1217      567591 :   const ulong ROOT8 = 251;
    1218      567591 :   const ulong ROOT7 = 563;
    1219      567591 :   const ulong ROOT5 = 7129;
    1220      567591 :   const ulong ROOT4 = 65521;
    1221             : #else
    1222       74585 :   const ulong ROOT9 = 11;
    1223       74585 :   const ulong ROOT8 = 13;
    1224       74585 :   const ulong ROOT7 = 23;
    1225       74585 :   const ulong ROOT5 = 83;
    1226       74585 :   const ulong ROOT4 = 251;
    1227             : #endif
    1228             :   ulong mask;
    1229             :   long v, i;
    1230             :   int e;
    1231      642176 :   if (n < 2) return 0;
    1232      642141 :   if (!odd(n)) {
    1233      341415 :     if (n & (n-1)) return 0;
    1234       65150 :     *pp = 2; return vals(n);
    1235             :   }
    1236      300726 :   if (n < 8) { *pp = n; return 1; } /* 3,5,7 */
    1237     3655068 :   for (i = 1/*skip p=2*/; i < TINYCUTOFF; i++)
    1238             :   {
    1239     3595993 :     ulong p = tinyprimes[i];
    1240     3595993 :     if (n % p == 0)
    1241             :     {
    1242      212426 :       v = u_lvalrem(n, p, &n);
    1243      212426 :       if (n == 1) { *pp = p; return v; }
    1244      209998 :       return 0;
    1245             :     }
    1246             :   }
    1247             :   /* p | n => p >= CUTOFF */
    1248             : 
    1249       59075 :   if (n < CUTOFF3)
    1250             :   {
    1251       46354 :     if (n < CUTOFF*CUTOFF || uisprime_101(n)) { *pp = n; return 1; }
    1252           0 :     if (uissquareall(n, &n)) { *pp = n; return 2; }
    1253           0 :     return 0;
    1254             :   }
    1255             : 
    1256             :   /* Check for squares, fourth powers, and eighth powers as appropriate. */
    1257       12721 :   v = 1;
    1258       12721 :   if (uissquareall(n, &n)) {
    1259           0 :     v <<= 1;
    1260           0 :     if (CUTOFF <= ROOT4 && uissquareall(n, &n)) {
    1261           0 :       v <<= 1;
    1262           0 :       if (CUTOFF <= ROOT8 && uissquareall(n, &n)) v <<= 1;
    1263             :     }
    1264             :   }
    1265             : 
    1266       12721 :   if (CUTOFF > ROOT5) mask = 1;
    1267             :   else
    1268             :   {
    1269       12720 :     const ulong CUTOFF5 = CUTOFF3*CUTOFF*CUTOFF;
    1270       12720 :     if (n < CUTOFF5) mask = 1; else mask = 3;
    1271       12720 :     if (CUTOFF <= ROOT7)
    1272             :     {
    1273       12720 :       const ulong CUTOFF7 = CUTOFF5*CUTOFF*CUTOFF;
    1274       12720 :       if (n >= CUTOFF7) mask = 7;
    1275             :     }
    1276             :   }
    1277             : 
    1278       12721 :   if (CUTOFF <= ROOT9 && (e = uis_357_power(n, &n, &mask))) { v *= e; mask=1; }
    1279       12721 :   if ((e = uis_357_power(n, &n, &mask))) v *= e;
    1280             : 
    1281       12721 :   if (uisprime_101(n)) { *pp = n; return v; }
    1282        6984 :   return 0;
    1283             : }
    1284             : 

Generated by: LCOV version 1.16