"Fossies" - the Fresh Open Source Software Archive

Member "pari-2.13.1/src/basemath/arith1.c" (4 Jan 2021, 154379 Bytes) of package /linux/misc/pari-2.13.1.tar.gz:


As a special service "Fossies" has tried to format the requested source page into HTML format using (guessed) C and C++ source code syntax highlighting (style: standard) with prefixed line numbers and code folding option. Alternatively you can here view or download the uninterpreted source code file. For more information about "arith1.c" see the Fossies "Dox" file reference documentation and the latest Fossies "Diffs" side-by-side code changes report: 2.13.0_vs_2.13.1.

    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 /*********************************************************************/
   16 /**                                                                 **/
   17 /**                     ARITHMETIC FUNCTIONS                        **/
   18 /**                         (first part)                            **/
   19 /**                                                                 **/
   20 /*********************************************************************/
   21 #include "pari.h"
   22 #include "paripriv.h"
   23 
   24 /******************************************************************/
   25 /*                                                                */
   26 /*                 GENERATOR of (Z/mZ)*                           */
   27 /*                                                                */
   28 /******************************************************************/
   29 static GEN
   30 remove2(GEN q) { long v = vali(q); return v? shifti(q, -v): q; }
   31 static ulong
   32 u_remove2(ulong q) { return q >> vals(q); }
   33 GEN
   34 odd_prime_divisors(GEN q) { return gel(Z_factor(remove2(q)), 1); }
   35 static GEN
   36 u_odd_prime_divisors(ulong q) { return gel(factoru(u_remove2(q)), 1); }
   37 /* p odd prime, q=(p-1)/2; L0 list of (some) divisors of q = (p-1)/2 or NULL
   38  * (all prime divisors of q); return the q/l, l in L0 */
   39 static GEN
   40 is_gener_expo(GEN p, GEN L0)
   41 {
   42   GEN L, q = shifti(p,-1);
   43   long i, l;
   44   if (L0) {
   45     l = lg(L0);
   46     L = cgetg(l, t_VEC);
   47   } else {
   48     L0 = L = odd_prime_divisors(q);
   49     l = lg(L);
   50   }
   51   for (i=1; i<l; i++) gel(L,i) = diviiexact(q, gel(L0,i));
   52   return L;
   53 }
   54 static GEN
   55 u_is_gener_expo(ulong p, GEN L0)
   56 {
   57   const ulong q = p >> 1;
   58   long i;
   59   GEN L;
   60   if (!L0) L0 = u_odd_prime_divisors(q);
   61   L = cgetg_copy(L0,&i);
   62   while (--i) L[i] = q / uel(L0,i);
   63   return L;
   64 }
   65 
   66 int
   67 is_gener_Fl(ulong x, ulong p, ulong p_1, GEN L)
   68 {
   69   long i;
   70   if (krouu(x, p) >= 0) return 0;
   71   for (i=lg(L)-1; i; i--)
   72   {
   73     ulong t = Fl_powu(x, uel(L,i), p);
   74     if (t == p_1 || t == 1) return 0;
   75   }
   76   return 1;
   77 }
   78 /* assume p prime */
   79 ulong
   80 pgener_Fl_local(ulong p, GEN L0)
   81 {
   82   const pari_sp av = avma;
   83   const ulong p_1 = p-1;
   84   long x;
   85   GEN L;
   86   if (p <= 19) switch(p)
   87   { /* quick trivial cases */
   88     case 2:  return 1;
   89     case 7:
   90     case 17: return 3;
   91     default: return 2;
   92   }
   93   L = u_is_gener_expo(p,L0);
   94   for (x = 2;; x++)
   95     if (is_gener_Fl(x,p,p_1,L)) return gc_ulong(av, x);
   96 }
   97 ulong
   98 pgener_Fl(ulong p) { return pgener_Fl_local(p, NULL); }
   99 
  100 /* L[i] = set of (p-1)/2l, l ODD prime divisor of p-1 (l=2 can be included,
  101  * but wasteful) */
  102 int
  103 is_gener_Fp(GEN x, GEN p, GEN p_1, GEN L)
  104 {
  105   long i, t = lgefint(x)==3? kroui(x[2], p): kronecker(x, p);
  106   if (t >= 0) return 0;
  107   for (i = lg(L)-1; i; i--)
  108   {
  109     GEN t = Fp_pow(x, gel(L,i), p);
  110     if (equalii(t, p_1) || equali1(t)) return 0;
  111   }
  112   return 1;
  113 }
  114 
  115 /* assume p prime, return a generator of all L[i]-Sylows in F_p^*. */
  116 GEN
  117 pgener_Fp_local(GEN p, GEN L0)
  118 {
  119   pari_sp av0 = avma;
  120   GEN x, p_1, L;
  121   if (lgefint(p) == 3)
  122   {
  123     ulong z;
  124     if (p[2] == 2) return gen_1;
  125     if (L0) L0 = ZV_to_nv(L0);
  126     z = pgener_Fl_local(uel(p,2), L0);
  127     set_avma(av0); return utoipos(z);
  128   }
  129   p_1 = subiu(p,1); L = is_gener_expo(p, L0);
  130   x = utoipos(2);
  131   for (;; x[2]++) { if (is_gener_Fp(x, p, p_1, L)) break; }
  132   set_avma(av0); return utoipos(uel(x,2));
  133 }
  134 
  135 GEN
  136 pgener_Fp(GEN p) { return pgener_Fp_local(p, NULL); }
  137 
  138 ulong
  139 pgener_Zl(ulong p)
  140 {
  141   if (p == 2) pari_err_DOMAIN("pgener_Zl","p","=",gen_2,gen_2);
  142   /* only p < 2^32 such that znprimroot(p) != znprimroot(p^2) */
  143   if (p == 40487) return 10;
  144 #ifndef LONG_IS_64BIT
  145   return pgener_Fl(p);
  146 #else
  147   if (p < (1UL<<32)) return pgener_Fl(p);
  148   else
  149   {
  150     const pari_sp av = avma;
  151     const ulong p_1 = p-1;
  152     long x ;
  153     GEN p2 = sqru(p), L = u_is_gener_expo(p, NULL);
  154     for (x=2;;x++)
  155       if (is_gener_Fl(x,p,p_1,L) && !is_pm1(Fp_powu(utoipos(x),p_1,p2)))
  156         return gc_ulong(av, x);
  157   }
  158 #endif
  159 }
  160 
  161 /* p prime. Return a primitive root modulo p^e, e > 1 */
  162 GEN
  163 pgener_Zp(GEN p)
  164 {
  165   if (lgefint(p) == 3) return utoipos(pgener_Zl(p[2]));
  166   else
  167   {
  168     const pari_sp av = avma;
  169     GEN p_1 = subiu(p,1), p2 = sqri(p), L = is_gener_expo(p,NULL);
  170     GEN x = utoipos(2);
  171     for (;; x[2]++)
  172       if (is_gener_Fp(x,p,p_1,L) && !equali1(Fp_pow(x,p_1,p2))) break;
  173     set_avma(av); return utoipos(uel(x,2));
  174   }
  175 }
  176 
  177 static GEN
  178 gener_Zp(GEN q, GEN F)
  179 {
  180   GEN p = NULL;
  181   long e = 0;
  182   if (F)
  183   {
  184     GEN P = gel(F,1), E = gel(F,2);
  185     long i, l = lg(P);
  186     for (i = 1; i < l; i++)
  187     {
  188       p = gel(P,i);
  189       if (absequaliu(p, 2)) continue;
  190       if (i < l-1) pari_err_DOMAIN("znprimroot", "argument","=",F,F);
  191       e = itos(gel(E,i));
  192     }
  193     if (!p) pari_err_DOMAIN("znprimroot", "argument","=",F,F);
  194   }
  195   else
  196     e = Z_isanypower(q, &p);
  197   return e > 1? pgener_Zp(p): pgener_Fp(q);
  198 }
  199 
  200 GEN
  201 znprimroot(GEN N)
  202 {
  203   pari_sp av = avma;
  204   GEN x, n, F;
  205 
  206   if ((F = check_arith_non0(N,"znprimroot")))
  207   {
  208     F = clean_Z_factor(F);
  209     N = typ(N) == t_VEC? gel(N,1): factorback(F);
  210   }
  211   N = absi_shallow(N);
  212   if (abscmpiu(N, 4) <= 0) { set_avma(av); return mkintmodu(N[2]-1,N[2]); }
  213   switch(mod4(N))
  214   {
  215     case 0: /* N = 0 mod 4 */
  216       pari_err_DOMAIN("znprimroot", "argument","=",N,N);
  217       x = NULL; break;
  218     case 2: /* N = 2 mod 4 */
  219       n = shifti(N,-1); /* becomes odd */
  220       x = gener_Zp(n,F); if (!mod2(x)) x = addii(x,n);
  221       break;
  222     default: /* N odd */
  223       x = gener_Zp(N,F);
  224       break;
  225   }
  226   return gerepilecopy(av, mkintmod(x, N));
  227 }
  228 
  229 /* n | (p-1), returns a primitive n-th root of 1 in F_p^* */
  230 GEN
  231 rootsof1_Fp(GEN n, GEN p)
  232 {
  233   pari_sp av = avma;
  234   GEN L = odd_prime_divisors(n); /* 2 implicit in pgener_Fp_local */
  235   GEN z = pgener_Fp_local(p, L);
  236   z = Fp_pow(z, diviiexact(subiu(p,1), n), p); /* prim. n-th root of 1 */
  237   return gerepileuptoint(av, z);
  238 }
  239 
  240 GEN
  241 rootsof1u_Fp(ulong n, GEN p)
  242 {
  243   pari_sp av = avma;
  244   GEN z, L = u_odd_prime_divisors(n); /* 2 implicit in pgener_Fp_local */
  245   z = pgener_Fp_local(p, Flv_to_ZV(L));
  246   z = Fp_pow(z, diviuexact(subiu(p,1), n), p); /* prim. n-th root of 1 */
  247   return gerepileuptoint(av, z);
  248 }
  249 
  250 ulong
  251 rootsof1_Fl(ulong n, ulong p)
  252 {
  253   pari_sp av = avma;
  254   GEN L = u_odd_prime_divisors(n); /* 2 implicit in pgener_Fl_local */
  255   ulong z = pgener_Fl_local(p, L);
  256   z = Fl_powu(z, (p-1) / n, p); /* prim. n-th root of 1 */
  257   return gc_ulong(av,z);
  258 }
  259 
  260 /*********************************************************************/
  261 /**                                                                 **/
  262 /**                     INVERSE TOTIENT FUNCTION                    **/
  263 /**                                                                 **/
  264 /*********************************************************************/
  265 /* N t_INT, L a ZV containing all prime divisors of N, and possibly other
  266  * primes. Return factor(N) */
  267 GEN
  268 Z_factor_listP(GEN N, GEN L)
  269 {
  270   long i, k, l = lg(L);
  271   GEN P = cgetg(l, t_COL), E = cgetg(l, t_COL);
  272   for (i = k = 1; i < l; i++)
  273   {
  274     GEN p = gel(L,i);
  275     long v = Z_pvalrem(N, p, &N);
  276     if (v)
  277     {
  278       gel(P,k) = p;
  279       gel(E,k) = utoipos(v);
  280       k++;
  281     }
  282   }
  283   setlg(P, k);
  284   setlg(E, k); return mkmat2(P,E);
  285 }
  286 
  287 /* look for x such that phi(x) = n, p | x => p > m (if m = NULL: no condition).
  288  * L is a list of primes containing all prime divisors of n. */
  289 static long
  290 istotient_i(GEN n, GEN m, GEN L, GEN *px)
  291 {
  292   pari_sp av = avma, av2;
  293   GEN k, D;
  294   long i, v;
  295   if (m && mod2(n))
  296   {
  297     if (!equali1(n)) return 0;
  298     if (px) *px = gen_1;
  299     return 1;
  300   }
  301   D = divisors(Z_factor_listP(shifti(n, -1), L));
  302   /* loop through primes p > m, d = p-1 | n */
  303   av2 = avma;
  304   if (!m)
  305   { /* special case p = 2, d = 1 */
  306     k = n;
  307     for (v = 1;; v++) {
  308       if (istotient_i(k, gen_2, L, px)) {
  309         if (px) *px = shifti(*px, v);
  310         return 1;
  311       }
  312       if (mod2(k)) break;
  313       k = shifti(k,-1);
  314     }
  315     set_avma(av2);
  316   }
  317   for (i = 1; i < lg(D); ++i)
  318   {
  319     GEN p, d = shifti(gel(D, i), 1); /* even divisors of n */
  320     if (m && cmpii(d, m) < 0) continue;
  321     p = addiu(d, 1);
  322     if (!isprime(p)) continue;
  323     k = diviiexact(n, d);
  324     for (v = 1;; v++) {
  325       GEN r;
  326       if (istotient_i(k, p, L, px)) {
  327         if (px) *px = mulii(*px, powiu(p, v));
  328         return 1;
  329       }
  330       k = dvmdii(k, p, &r);
  331       if (r != gen_0) break;
  332     }
  333     set_avma(av2);
  334   }
  335   return gc_long(av,0);
  336 }
  337 
  338 /* find x such that phi(x) = n */
  339 long
  340 istotient(GEN n, GEN *px)
  341 {
  342   pari_sp av = avma;
  343   if (typ(n) != t_INT) pari_err_TYPE("istotient", n);
  344   if (signe(n) < 1) return 0;
  345   if (mod2(n))
  346   {
  347     if (!equali1(n)) return 0;
  348     if (px) *px = gen_1;
  349     return 1;
  350   }
  351   if (istotient_i(n, NULL, gel(Z_factor(n), 1), px))
  352   {
  353     if (!px) set_avma(av);
  354     else
  355       *px = gerepileuptoint(av, *px);
  356     return 1;
  357   }
  358   return gc_long(av,0);
  359 }
  360 
  361 /*********************************************************************/
  362 /**                                                                 **/
  363 /**                     INTEGRAL LOGARITHM                          **/
  364 /**                                                                 **/
  365 /*********************************************************************/
  366 
  367 /* y > 1 and B > 0 integers. Return e such that y^e <= B < y^(e+1), i.e
  368  * e = floor(log_y B). Set *ptq = y^e if non-NULL */
  369 long
  370 ulogintall(ulong B, ulong y, ulong *ptq)
  371 {
  372   ulong r, r2;
  373   long e;
  374 
  375   if (y == 2)
  376   {
  377     long eB = expu(B); /* 2^eB <= B < 2^(eB + 1) */
  378     if (ptq) *ptq = 1UL << eB;
  379     return eB;
  380   }
  381   r = y, r2 = 1UL;
  382   for (e=1;; e++)
  383   { /* here, r = y^e, r2 = y^(e-1) */
  384     if (r >= B)
  385     {
  386       if (r != B) { e--; r = r2; }
  387       if (ptq) *ptq = r;
  388       return e;
  389     }
  390     r2 = r;
  391     r = umuluu_or_0(y, r);
  392     if (!r)
  393     {
  394       if (ptq) *ptq = r2;
  395       return e;
  396     }
  397   }
  398 }
  399 
  400 /* y > 1 and B > 0 integers. Return e such that y^e <= B < y^(e+1), i.e
  401  * e = floor(log_y B). Set *ptq = y^e if non-NULL */
  402 long
  403 logintall(GEN B, GEN y, GEN *ptq)
  404 {
  405   pari_sp av;
  406   long ey, e, emax, i, eB = expi(B); /* 2^eB <= B < 2^(eB + 1) */
  407   GEN q, pow2;
  408 
  409   if (lgefint(B) == 3)
  410   {
  411     ulong q;
  412     if (lgefint(y) > 3)
  413     {
  414       if (ptq) *ptq = gen_1;
  415       return 0;
  416     }
  417     if (!ptq) return ulogintall(B[2], y[2], NULL);
  418     e = ulogintall(B[2], y[2], &q);
  419     *ptq = utoi(q); return e;
  420   }
  421   if (equaliu(y,2))
  422   {
  423     if (ptq) *ptq = int2n(eB);
  424     return eB;
  425   }
  426   av = avma;
  427   ey = expi(y);
  428   /* eB/(ey+1) - 1 < e <= eB/ey */
  429   emax = eB/ey;
  430   if (emax <= 13) /* e small, be naive */
  431   {
  432     GEN r = y, r2 = gen_1;
  433     for (e=1;; e++)
  434     { /* here, r = y^e, r2 = y^(e-1) */
  435       long fl = cmpii(r, B);
  436       if (fl >= 0)
  437       {
  438         if (fl) { e--; cgiv(r); r = r2; }
  439         if (ptq) *ptq = gerepileuptoint(av, r); else set_avma(av);
  440         return e;
  441       }
  442       r2 = r; r = mulii(r,y);
  443     }
  444   }
  445   /* e >= 13 ey / (ey+1) >= 6.5 */
  446 
  447   /* binary splitting: compute bits of e one by one */
  448   /* compute pow2[i] = y^(2^i) [i < crude upper bound for log_2 log_y(B)] */
  449   pow2 = new_chunk((long)log2(eB)+2);
  450   gel(pow2,0) = y;
  451   for (i=0, q=y;; )
  452   {
  453     GEN r = gel(pow2,i); /* r = y^2^i */
  454     long fl = cmpii(r,B);
  455     if (!fl)
  456     {
  457       e = 1L<<i;
  458       if (ptq) *ptq = gerepileuptoint(av, r); else set_avma(av);
  459       return e;
  460     }
  461     if (fl > 0) { i--; break; }
  462     q = r;
  463     if (1L<<(i+1) > emax) break;
  464     gel(pow2,++i) = sqri(q);
  465   }
  466 
  467   for (e = 1L<<i;;)
  468   { /* y^e = q < B < r = q * y^(2^i) */
  469     pari_sp av2 = avma;
  470     long fl;
  471     GEN r;
  472     if (--i < 0) break;
  473     r = mulii(q, gel(pow2,i));
  474     fl = cmpii(r, B);
  475     if (fl > 0) set_avma(av2);
  476     else
  477     {
  478       e += (1L<<i);
  479       q = r;
  480       if (!fl) break; /* B = r */
  481     }
  482   }
  483   if (ptq) *ptq = gerepileuptoint(av, q); else set_avma(av);
  484   return e;
  485 }
  486 
  487 long
  488 logint0(GEN B, GEN y, GEN *ptq)
  489 {
  490   if (typ(B) != t_INT) pari_err_TYPE("logint",B);
  491   if (signe(B) <= 0) pari_err_DOMAIN("logint", "x" ,"<=", gen_0, B);
  492   if (typ(y) != t_INT) pari_err_TYPE("logint",y);
  493   if (cmpis(y, 2) < 0) pari_err_DOMAIN("logint", "b" ,"<=", gen_1, y);
  494   return logintall(B,y,ptq);
  495 }
  496 
  497 /*********************************************************************/
  498 /**                                                                 **/
  499 /**                     INTEGRAL SQUARE ROOT                        **/
  500 /**                                                                 **/
  501 /*********************************************************************/
  502 GEN
  503 sqrtint(GEN a)
  504 {
  505   if (typ(a) != t_INT) pari_err_TYPE("sqrtint",a);
  506   switch (signe(a))
  507   {
  508     case 1: return sqrti(a);
  509     case 0: return gen_0;
  510     default: pari_err_DOMAIN("sqrtint", "argument", "<", gen_0,a);
  511   }
  512   return NULL; /* LCOV_EXCL_LINE */
  513 }
  514 GEN
  515 sqrtint0(GEN a, GEN *r)
  516 {
  517   if (!r) return sqrtint(a);
  518   if (typ(a) != t_INT) pari_err_TYPE("sqrtint",a);
  519   switch (signe(a))
  520   {
  521     case 1: return sqrtremi(a, r);
  522     case 0: *r = gen_0; return gen_0;
  523     default: pari_err_DOMAIN("sqrtint", "argument", "<", gen_0,a);
  524   }
  525   return NULL; /* LCOV_EXCL_LINE */
  526 }
  527 
  528 /*********************************************************************/
  529 /**                                                                 **/
  530 /**                      PERFECT SQUARE                             **/
  531 /**                                                                 **/
  532 /*********************************************************************/
  533 static int
  534 carremod(ulong A)
  535 {
  536   const int carresmod64[]={
  537     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,
  538     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};
  539   const int carresmod63[]={
  540     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,
  541     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};
  542   const int carresmod65[]={
  543     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,
  544     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};
  545   const int carresmod11[]={1,1,0,1,1,1,0,0,0,1, 0};
  546   return (carresmod64[A & 0x3fUL]
  547     && carresmod63[A % 63UL]
  548     && carresmod65[A % 65UL]
  549     && carresmod11[A % 11UL]);
  550 }
  551 
  552 /* emulate Z_issquareall on single-word integers */
  553 long
  554 uissquareall(ulong A, ulong *sqrtA)
  555 {
  556   if (!A) { *sqrtA = 0; return 1; }
  557   if (carremod(A))
  558   {
  559     ulong a = usqrt(A);
  560     if (a * a == A) { *sqrtA = a; return 1; }
  561   }
  562   return 0;
  563 }
  564 long
  565 uissquare(ulong A)
  566 {
  567   if (!A) return 1;
  568   if (carremod(A))
  569   {
  570     ulong a = usqrt(A);
  571     if (a * a == A) return 1;
  572   }
  573   return 0;
  574 }
  575 
  576 long
  577 Z_issquareall(GEN x, GEN *pt)
  578 {
  579   pari_sp av;
  580   GEN y, r;
  581 
  582   switch(signe(x))
  583   {
  584     case -1: return 0;
  585     case 0: if (pt) *pt=gen_0; return 1;
  586   }
  587   if (lgefint(x) == 3)
  588   {
  589     ulong u = uel(x,2), a;
  590     if (!pt) return uissquare(u);
  591     if (!uissquareall(u, &a)) return 0;
  592     *pt = utoipos(a); return 1;
  593   }
  594   if (!carremod(umodiu(x, 64*63*65*11))) return 0;
  595   av = avma; y = sqrtremi(x, &r);
  596   if (r != gen_0) return gc_long(av,0);
  597   if (pt) { *pt = y; set_avma((pari_sp)y); } else set_avma(av);
  598   return 1;
  599 }
  600 
  601 /* a t_INT, p prime */
  602 long
  603 Zp_issquare(GEN a, GEN p)
  604 {
  605   long v;
  606   GEN ap;
  607 
  608   if (!signe(a) || gequal1(a)) return 1;
  609   v = Z_pvalrem(a, p, &ap);
  610   if (v&1) return 0;
  611   return absequaliu(p, 2)? umodiu(ap, 8) == 1
  612                       : kronecker(ap,p) == 1;
  613 }
  614 
  615 static long
  616 polissquareall(GEN x, GEN *pt)
  617 {
  618   pari_sp av;
  619   long v;
  620   GEN y, a, b, p;
  621 
  622   if (!signe(x))
  623   {
  624     if (pt) *pt = gcopy(x);
  625     return 1;
  626   }
  627   if (odd(degpol(x))) return 0; /* odd degree */
  628   av = avma;
  629   v = RgX_valrem(x, &x);
  630   if (v & 1) return gc_long(av,0);
  631   a = gel(x,2); /* test constant coeff */
  632   if (!pt)
  633   { if (!issquare(a)) return gc_long(av,0); }
  634   else
  635   { if (!issquareall(a,&b)) return gc_long(av,0); }
  636   if (!degpol(x)) { /* constant polynomial */
  637     if (!pt) return gc_long(av,1);
  638     y = scalarpol(b, varn(x)); goto END;
  639   }
  640   p = characteristic(x);
  641   if (signe(p) && !mod2(p))
  642   {
  643     long i, lx;
  644     if (!absequaliu(p,2)) pari_err_IMPL("issquare for even characteristic != 2");
  645     x = gmul(x, mkintmod(gen_1, gen_2));
  646     lx = lg(x);
  647     if ((lx-3) & 1) return gc_long(av,0);
  648     for (i = 3; i < lx; i+=2)
  649       if (!gequal0(gel(x,i))) return gc_long(av,0);
  650     if (pt) {
  651       y = cgetg((lx+3) / 2, t_POL);
  652       for (i = 2; i < lx; i+=2)
  653         if (!issquareall(gel(x,i), &gel(y, (i+2)>>1))) return gc_long(av,0);
  654       y[1] = evalsigne(1) | evalvarn(varn(x));
  655       goto END;
  656     } else {
  657       for (i = 2; i < lx; i+=2)
  658         if (!issquare(gel(x,i))) return gc_long(av,0);
  659       return gc_long(av,1);
  660     }
  661   }
  662   else
  663   {
  664     long m = 1;
  665     x = RgX_Rg_div(x,a);
  666     /* a(x^m) = B^2 => B = b(x^m) provided a(0) != 0 */
  667     if (!signe(p)) x = RgX_deflate_max(x,&m);
  668     y = ser2rfrac_i(gsqrt(RgX_to_ser(x,lg(x)-1),0));
  669     if (!RgX_equal(RgX_sqr(y), x)) return gc_long(av,0);
  670     if (!pt) return gc_long(av,1);
  671     if (!gequal1(a)) y = gmul(b, y);
  672     if (m != 1) y = RgX_inflate(y,m);
  673   }
  674 END:
  675   if (v) y = RgX_shift_shallow(y, v>>1);
  676   *pt = gerepilecopy(av, y); return 1;
  677 }
  678 
  679 /* b unit mod p */
  680 static int
  681 Up_ispower(GEN b, GEN K, GEN p, long d, GEN *pt)
  682 {
  683   if (d == 1)
  684   { /* mod p: faster */
  685     if (!Fp_ispower(b, K, p)) return 0;
  686     if (pt) *pt = Fp_sqrtn(b, K, p, NULL);
  687   }
  688   else
  689   { /* mod p^{2 +} */
  690     if (!ispower(cvtop(b, p, d), K, pt)) return 0;
  691     if (pt) *pt = gtrunc(*pt);
  692   }
  693   return 1;
  694 }
  695 
  696 /* We're studying whether a mod (q*p^e) is a K-th power, (q,p) = 1.
  697  * Decide mod p^e, then reduce a mod q unless q = NULL. */
  698 static int
  699 handle_pe(GEN *pa, GEN q, GEN L, GEN K, GEN p, long e)
  700 {
  701   GEN t, A;
  702   long v = Z_pvalrem(*pa, p, &A), d = e - v;
  703   if (d <= 0) t = gen_0;
  704   else
  705   {
  706     ulong r;
  707     v = uabsdivui_rem(v, K, &r);
  708     if (r || !Up_ispower(A, K, p, d, L? &t: NULL)) return 0;
  709     if (L && v) t = mulii(t, powiu(p, v));
  710   }
  711   if (q) *pa = modii(*pa, q);
  712   if (L) vectrunc_append(L, mkintmod(t, powiu(p, e)));
  713   return 1;
  714 }
  715 long
  716 Zn_ispower(GEN a, GEN q, GEN K, GEN *pt)
  717 {
  718   GEN L, N;
  719   pari_sp av;
  720   long e, i, l;
  721   ulong pp;
  722   forprime_t S;
  723 
  724   if (!signe(a))
  725   {
  726     if (pt) {
  727       GEN t = cgetg(3, t_INTMOD);
  728       gel(t,1) = icopy(q); gel(t,2) = gen_0; *pt = t;
  729     }
  730     return 1;
  731   }
  732   /* a != 0 */
  733   av = avma;
  734 
  735   if (typ(q) != t_INT) /* integer factorization */
  736   {
  737     GEN P = gel(q,1), E = gel(q,2);
  738     l = lg(P);
  739     L = pt? vectrunc_init(l): NULL;
  740     for (i = 1; i < l; i++)
  741     {
  742       GEN p = gel(P,i);
  743       long e = itos(gel(E,i));
  744       if (!handle_pe(&a, NULL, L, K, p, e)) return gc_long(av,0);
  745     }
  746     goto END;
  747   }
  748   if (!mod2(K)
  749       && kronecker(a, shifti(q,-vali(q))) == -1) return gc_long(av,0);
  750   L = pt? vectrunc_init(expi(q)+1): NULL;
  751   u_forprime_init(&S, 2, tridiv_bound(q));
  752   while ((pp = u_forprime_next(&S)))
  753   {
  754     int stop;
  755     e = Z_lvalrem_stop(&q, pp, &stop);
  756     if (!e) continue;
  757     if (!handle_pe(&a, q, L, K, utoipos(pp), e)) return gc_long(av,0);
  758     if (stop)
  759     {
  760       if (!is_pm1(q) && !handle_pe(&a, q, L, K, q, 1)) return gc_long(av,0);
  761       goto END;
  762     }
  763   }
  764   l = lg(primetab);
  765   for (i = 1; i < l; i++)
  766   {
  767     GEN p = gel(primetab,i);
  768     e = Z_pvalrem(q, p, &q);
  769     if (!e) continue;
  770     if (!handle_pe(&a, q, L, K, p, e)) return gc_long(av,0);
  771     if (is_pm1(q)) goto END;
  772   }
  773   N = gcdii(a,q);
  774   if (!is_pm1(N))
  775   {
  776     if (ifac_isprime(N))
  777     {
  778       e = Z_pvalrem(q, N, &q);
  779       if (!handle_pe(&a, q, L, K, N, e)) return gc_long(av,0);
  780     }
  781     else
  782     {
  783       GEN part = ifac_start(N, 0);
  784       for(;;)
  785       {
  786         long e;
  787         GEN p;
  788         if (!ifac_next(&part, &p, &e)) break;
  789         e = Z_pvalrem(q, p, &q);
  790         if (!handle_pe(&a, q, L, K, p, e)) return gc_long(av,0);
  791       }
  792     }
  793   }
  794   if (!is_pm1(q))
  795   {
  796     if (ifac_isprime(q))
  797     {
  798       if (!handle_pe(&a, q, L, K, q, 1)) return gc_long(av,0);
  799     }
  800     else
  801     {
  802       GEN part = ifac_start(q, 0);
  803       for(;;)
  804       {
  805         long e;
  806         GEN p;
  807         if (!ifac_next(&part, &p, &e)) break;
  808         if (!handle_pe(&a, q, L, K, p, e)) return gc_long(av,0);
  809       }
  810     }
  811   }
  812 END:
  813   if (pt) *pt = gerepileupto(av, chinese1_coprime_Z(L));
  814   return 1;
  815 }
  816 
  817 static long
  818 polmodispower(GEN x, GEN K, GEN *pt)
  819 {
  820   pari_sp av = avma;
  821   GEN p = NULL, T = NULL;
  822   if (Rg_is_FpXQ(x, &T,&p) && p)
  823   {
  824     x = liftall_shallow(x);
  825     if (T) T = liftall_shallow(T);
  826     if (!Fq_ispower(x, K, T, p)) return gc_long(av,0);
  827     if (!pt) return gc_long(av,1);
  828     x = Fq_sqrtn(x, K, T,p, NULL);
  829     if (typ(x) == t_INT)
  830       x = Fp_to_mod(x,p);
  831     else
  832       x = mkpolmod(FpX_to_mod(x,p), FpX_to_mod(T,p));
  833     *pt = gerepilecopy(av, x); return 1;
  834   }
  835   pari_err_IMPL("ispower for general t_POLMOD");
  836   return 0;
  837 }
  838 
  839 long
  840 issquareall(GEN x, GEN *pt)
  841 {
  842   long tx = typ(x);
  843   GEN F;
  844   pari_sp av;
  845 
  846   if (!pt) return issquare(x);
  847   switch(tx)
  848   {
  849     case t_INT: return Z_issquareall(x, pt);
  850     case t_FRAC: av = avma;
  851       F = cgetg(3, t_FRAC);
  852       if (   !Z_issquareall(gel(x,1), &gel(F,1))
  853           || !Z_issquareall(gel(x,2), &gel(F,2))) return gc_long(av,0);
  854       *pt = F; return 1;
  855 
  856     case t_POLMOD:
  857       return polmodispower(x, gen_2, pt);
  858     case t_POL: return polissquareall(x,pt);
  859     case t_RFRAC: av = avma;
  860       F = cgetg(3, t_RFRAC);
  861       if (   !issquareall(gel(x,1), &gel(F,1))
  862           || !polissquareall(gel(x,2), &gel(F,2))) return gc_long(av,0);
  863       *pt = F; return 1;
  864 
  865     case t_REAL: case t_COMPLEX: case t_PADIC: case t_SER:
  866       if (!issquare(x)) return 0;
  867       *pt = gsqrt(x, DEFAULTPREC); return 1;
  868 
  869     case t_INTMOD:
  870       return Zn_ispower(gel(x,2), gel(x,1), gen_2, pt);
  871 
  872     case t_FFELT: return FF_issquareall(x, pt);
  873 
  874   }
  875   pari_err_TYPE("issquareall",x);
  876   return 0; /* LCOV_EXCL_LINE */
  877 }
  878 
  879 long
  880 issquare(GEN x)
  881 {
  882   pari_sp av;
  883   GEN a, p;
  884   long v;
  885 
  886   switch(typ(x))
  887   {
  888     case t_INT:
  889       return Z_issquare(x);
  890 
  891     case t_REAL:
  892       return (signe(x)>=0);
  893 
  894     case t_INTMOD:
  895       return Zn_ispower(gel(x,2), gel(x,1), gen_2, NULL);
  896 
  897     case t_FRAC:
  898       return Z_issquare(gel(x,1)) && Z_issquare(gel(x,2));
  899 
  900     case t_FFELT: return FF_issquareall(x, NULL);
  901 
  902     case t_COMPLEX:
  903       return 1;
  904 
  905     case t_PADIC:
  906       a = gel(x,4); if (!signe(a)) return 1;
  907       if (valp(x)&1) return 0;
  908       p = gel(x,2);
  909       if (!absequaliu(p, 2)) return (kronecker(a,p) != -1);
  910 
  911       v = precp(x); /* here p=2, a is odd */
  912       if ((v>=3 && mod8(a) != 1 ) ||
  913           (v==2 && mod4(a) != 1)) return 0;
  914       return 1;
  915 
  916     case t_POLMOD:
  917       return polmodispower(x, gen_2, NULL);
  918 
  919     case t_POL:
  920       return polissquareall(x,NULL);
  921 
  922     case t_SER:
  923       if (!signe(x)) return 1;
  924       if (valp(x)&1) return 0;
  925       return issquare(gel(x,2));
  926 
  927     case t_RFRAC:
  928       av = avma; return gc_long(av, issquare(gmul(gel(x,1),gel(x,2))));
  929   }
  930   pari_err_TYPE("issquare",x);
  931   return 0; /* LCOV_EXCL_LINE */
  932 }
  933 GEN gissquare(GEN x) { return issquare(x)? gen_1: gen_0; }
  934 GEN gissquareall(GEN x, GEN *pt) { return issquareall(x,pt)? gen_1: gen_0; }
  935 
  936 long
  937 ispolygonal(GEN x, GEN S, GEN *N)
  938 {
  939   pari_sp av = avma;
  940   GEN D, d, n;
  941   if (typ(x) != t_INT) pari_err_TYPE("ispolygonal", x);
  942   if (typ(S) != t_INT) pari_err_TYPE("ispolygonal", S);
  943   if (abscmpiu(S,3) < 0) pari_err_DOMAIN("ispolygonal","s","<", utoipos(3),S);
  944   if (signe(x) < 0) return 0;
  945   if (signe(x) == 0) { if (N) *N = gen_0; return 1; }
  946   if (is_pm1(x)) { if (N) *N = gen_1; return 1; }
  947   /* n = (sqrt( (8s - 16) x + (s-4)^2 ) + s - 4) / 2(s - 2) */
  948   if (abscmpiu(S, 1<<16) < 0) /* common case ! */
  949   {
  950     ulong s = S[2], r;
  951     if (s == 4) return Z_issquareall(x, N);
  952     if (s == 3)
  953       D = addiu(shifti(x, 3), 1);
  954     else
  955       D = addiu(mului(8*s - 16, x), (s-4)*(s-4));
  956     if (!Z_issquareall(D, &d)) return gc_long(av,0);
  957     if (s == 3)
  958       d = subiu(d, 1);
  959     else
  960       d = addiu(d, s - 4);
  961     n = absdiviu_rem(d, 2*s - 4, &r);
  962     if (r) return gc_long(av,0);
  963   }
  964   else
  965   {
  966     GEN r, S_2 = subiu(S,2), S_4 = subiu(S,4);
  967     D = addii(mulii(shifti(S_2,3), x), sqri(S_4));
  968     if (!Z_issquareall(D, &d)) return gc_long(av,0);
  969     d = addii(d, S_4);
  970     n = dvmdii(shifti(d,-1), S_2, &r);
  971     if (r != gen_0) return gc_long(av,0);
  972   }
  973   if (N) *N = gerepileuptoint(av, n); else set_avma(av);
  974   return 1;
  975 }
  976 
  977 /*********************************************************************/
  978 /**                                                                 **/
  979 /**                        PERFECT POWER                            **/
  980 /**                                                                 **/
  981 /*********************************************************************/
  982 static long
  983 polispower(GEN x, GEN K, GEN *pt)
  984 {
  985   pari_sp av;
  986   long v, d, k = itos(K);
  987   GEN y, a, b;
  988   GEN T = NULL, p = NULL;
  989 
  990   if (!signe(x))
  991   {
  992     if (pt) *pt = gcopy(x);
  993     return 1;
  994   }
  995   d = degpol(x);
  996   if (d % k) return 0; /* degree not multiple of k */
  997   av = avma;
  998   if (RgX_is_FpXQX(x, &T, &p) && p)
  999   { /* over Fq */
 1000     if (T && typ(T) == t_FFELT)
 1001     {
 1002       if (!FFX_ispower(x, k, T, pt)) return gc_long(av,0);
 1003       return 1;
 1004     }
 1005     x = RgX_to_FqX(x,T,p);
 1006     if (!FqX_ispower(x, k, T,p, pt)) return gc_long(av,0);
 1007     if (pt) *pt = gerepileupto(av, FqX_to_mod(*pt, T, p));
 1008     return 1;
 1009   }
 1010   v = RgX_valrem(x, &x);
 1011   if (v % k) return 0;
 1012   v /= k;
 1013   a = gel(x,2); b = NULL;
 1014   if (!ispower(a, K, &b)) return gc_long(av,0);
 1015   if (d)
 1016   {
 1017     GEN p = characteristic(x);
 1018     a = leading_coeff(x);
 1019     if (!ispower(a, K, &b)) return gc_long(av,0);
 1020     x = RgX_normalize(x);
 1021     if (signe(p) && cmpii(p,K) <= 0)
 1022       pari_err_IMPL("ispower(general t_POL) in small characteristic");
 1023     y = gtrunc(gsqrtn(RgX_to_ser(x,lg(x)), K, NULL, 0));
 1024     if (!RgX_equal(powgi(y, K), x)) return gc_long(av,0);
 1025   }
 1026   else
 1027     y = pol_1(varn(x));
 1028   if (pt)
 1029   {
 1030     if (!gequal1(a))
 1031     {
 1032       if (!b) b = gsqrtn(a, K, NULL, DEFAULTPREC);
 1033       y = gmul(b,y);
 1034     }
 1035     if (v) y = RgX_shift_shallow(y, v);
 1036     *pt = gerepilecopy(av, y);
 1037   }
 1038   else set_avma(av);
 1039   return 1;
 1040 }
 1041 
 1042 long
 1043 Z_ispowerall(GEN x, ulong k, GEN *pt)
 1044 {
 1045   long s = signe(x);
 1046   ulong mask;
 1047   if (!s) { if (pt) *pt = gen_0; return 1; }
 1048   if (s > 0) {
 1049     if (k == 2) return Z_issquareall(x, pt);
 1050     if (k == 3) { mask = 1; return !!is_357_power(x, pt, &mask); }
 1051     if (k == 5) { mask = 2; return !!is_357_power(x, pt, &mask); }
 1052     if (k == 7) { mask = 4; return !!is_357_power(x, pt, &mask); }
 1053     return is_kth_power(x, k, pt);
 1054   }
 1055   if (!odd(k)) return 0;
 1056   if (Z_ispowerall(absi_shallow(x), k, pt))
 1057   {
 1058     if (pt) *pt = negi(*pt);
 1059     return 1;
 1060   };
 1061   return 0;
 1062 }
 1063 
 1064 /* is x a K-th power mod p ? Assume p prime. */
 1065 int
 1066 Fp_ispower(GEN x, GEN K, GEN p)
 1067 {
 1068   pari_sp av = avma;
 1069   GEN p_1;
 1070   x = modii(x, p);
 1071   if (!signe(x) || equali1(x)) return gc_bool(av,1);
 1072   /* implies p > 2 */
 1073   p_1 = subiu(p,1);
 1074   K = gcdii(K, p_1);
 1075   if (absequaliu(K, 2)) return gc_bool(av, kronecker(x,p) > 0);
 1076   x = Fp_pow(x, diviiexact(p_1,K), p);
 1077   return gc_bool(av, equali1(x));
 1078 }
 1079 
 1080 /* x unit defined modulo 2^e, e > 0, p prime */
 1081 static int
 1082 U2_issquare(GEN x, long e)
 1083 {
 1084   long r = signe(x)>=0?mod8(x):8-mod8(x);
 1085   if (e==1) return 1;
 1086   if (e==2) return (r&3L) == 1;
 1087   return r == 1;
 1088 }
 1089 /* x unit defined modulo p^e, e > 0, p prime */
 1090 static int
 1091 Up_issquare(GEN x, GEN p, long e)
 1092 { return (absequaliu(p,2))? U2_issquare(x, e): kronecker(x,p)==1; }
 1093 
 1094 long
 1095 Zn_issquare(GEN d, GEN fn)
 1096 {
 1097   long j, np;
 1098   if (typ(d) != t_INT) pari_err_TYPE("Zn_issquare",d);
 1099   if (typ(fn) == t_INT) return Zn_ispower(d, fn, gen_2, NULL);
 1100   /* integer factorization */
 1101   np = nbrows(fn);
 1102   for (j = 1; j <= np; ++j)
 1103   {
 1104     GEN  r, p = gcoeff(fn, j, 1);
 1105     long e = itos(gcoeff(fn, j, 2));
 1106     long v = Z_pvalrem(d,p,&r);
 1107     if (v < e && (odd(v) || !Up_issquare(r, p, e-v))) return 0;
 1108   }
 1109   return 1;
 1110 }
 1111 
 1112 /* return [N',v]; v contains all x mod N' s.t. x^2 + B x + C = 0 modulo N */
 1113 GEN
 1114 Zn_quad_roots(GEN N, GEN B, GEN C)
 1115 {
 1116   pari_sp av = avma;
 1117   GEN fa = NULL, D, w, v, P, E, F0, Q0, F, mF, A, Q, T, R, Np, N4;
 1118   long l, i, j, ct;
 1119 
 1120   if ((fa = check_arith_non0(N,"Zn_quad_roots")))
 1121   {
 1122     N = typ(N) == t_VEC? gel(N,1): factorback(N);
 1123     fa = clean_Z_factor(fa);
 1124   }
 1125   N = absi_shallow(N);
 1126   N4 = shifti(N,2);
 1127   D = modii(subii(sqri(B), shifti(C,2)), N4);
 1128   if (!signe(D))
 1129   { /* (x + B/2)^2 = 0 (mod N), D = B^2-4C = 0 (4N) => B even */
 1130     if (!fa) fa = Z_factor(N);
 1131     P = gel(fa,1);
 1132     E = ZV_to_zv(gel(fa,2));
 1133     l = lg(P);
 1134     for (i = 1; i < l; i++) E[i] = (E[i]+1) >> 1;
 1135     Np = factorback2(P, E); /* x = -B mod N' */
 1136     B = shifti(B,-1);
 1137     return gerepilecopy(av, mkvec2(Np, mkvec(Fp_neg(B,Np))));
 1138   }
 1139   if (!fa)
 1140     fa = Z_factor(N4);
 1141   else  /* convert to factorization of N4 = 4*N */
 1142     fa = famat_reduce(famat_mulpows_shallow(fa, gen_2, 2));
 1143   P = gel(fa,1); l = lg(P);
 1144   E = ZV_to_zv(gel(fa,2));
 1145   F = cgetg(l, t_VEC);
 1146   mF= cgetg(l, t_VEC); F0 = gen_0;
 1147   Q = cgetg(l, t_VEC); Q0 = gen_1;
 1148   for (i = j = 1, ct = 0; i < l; i++)
 1149   {
 1150     GEN p = gel(P,i), q, f, mf, D0;
 1151     long t2, s = E[i], t = Z_pvalrem(D, p, &D0), d = s - t;
 1152     if (d <= 0)
 1153     {
 1154       q = powiu(p, (s+1)>>1);
 1155       Q0 = mulii(Q0, q); continue;
 1156     }
 1157     /* d > 0 */
 1158     if (odd(t)) return NULL;
 1159     t2 = t >> 1;
 1160     if (i > 1)
 1161     { /* p > 2 */
 1162       if (kronecker(D0, p) == -1) return NULL;
 1163       q = powiu(p, s - t2);
 1164       f = Zp_sqrt(D0, p, d);
 1165       if (!f) return NULL; /* p was not actually prime... */
 1166       if (t2) f = mulii(powiu(p,t2), f);
 1167       mf = Fp_neg(f, q);
 1168     }
 1169     else
 1170     { /* p = 2 */
 1171       if (d <= 3)
 1172       {
 1173         if (d == 3 && Mod8(D0) != 1) return NULL;
 1174         if (d == 2 && Mod4(D0) != 1) return NULL;
 1175         Q0 = int2n(1+t2); F0 = NULL; continue;
 1176       }
 1177       if (Mod8(D0) != 1) return NULL;
 1178       q = int2n(d - 1 + t2);
 1179       f = Z2_sqrt(D0, d);
 1180       if (t2) f = shifti(f, t2);
 1181       mf = Fp_neg(f, q);
 1182     }
 1183     gel(Q,j) = q;
 1184     gel(F,j) = f;
 1185     gel(mF,j)= mf; j++;
 1186   }
 1187   setlg(Q,j);
 1188   setlg(F,j);
 1189   setlg(mF,j);
 1190   if (is_pm1(Q0)) A = leafcopy(F);
 1191   else
 1192   { /* append the fixed congruence (F0 mod Q0) */
 1193     if (!F0) F0 = shifti(Q0,-1);
 1194     A = shallowconcat(F, F0);
 1195     Q = shallowconcat(Q, Q0);
 1196   }
 1197   ct = 1 << (j-1);
 1198   T = ZV_producttree(Q);
 1199   R = ZV_chinesetree(Q,T);
 1200   Np = gmael(T, lg(T)-1, 1);
 1201   B = modii(B, Np);
 1202   if (!signe(B)) B = NULL;
 1203   Np = shifti(Np, -1); /* N' = (\prod_i Q[i]) / 2 */
 1204   w = cgetg(3, t_VEC);
 1205   gel(w,1) = icopy(Np);
 1206   gel(w,2) = v = cgetg(ct+1, t_VEC);
 1207   l = lg(F);
 1208   for (j = 1; j <= ct; j++)
 1209   {
 1210     pari_sp av2 = avma;
 1211     long m = j - 1;
 1212     GEN u;
 1213     for (i = 1; i < l; i++)
 1214     {
 1215       gel(A,i) = (m&1L)? gel(mF,i): gel(F,i);
 1216       m >>= 1;
 1217     }
 1218     u = ZV_chinese_tree(A,Q,T,R); /* u mod N' st u^2 = B^2-4C modulo 4N */
 1219     if (B) u = subii(u,B);
 1220     gel(v,j) = gerepileuptoint(av2, modii(shifti(u,-1), Np));
 1221   }
 1222   return gerepileupto(av, w);
 1223 }
 1224 
 1225 static long
 1226 Qp_ispower(GEN x, GEN K, GEN *pt)
 1227 {
 1228   pari_sp av = avma;
 1229   GEN z = Qp_sqrtn(x, K, NULL);
 1230   if (!z) return gc_long(av,0);
 1231   if (pt) *pt = z;
 1232   return 1;
 1233 }
 1234 
 1235 long
 1236 ispower(GEN x, GEN K, GEN *pt)
 1237 {
 1238   GEN z;
 1239 
 1240   if (!K) return gisanypower(x, pt);
 1241   if (typ(K) != t_INT) pari_err_TYPE("ispower",K);
 1242   if (signe(K) <= 0) pari_err_DOMAIN("ispower","exponent","<=",gen_0,K);
 1243   if (equali1(K)) { if (pt) *pt = gcopy(x); return 1; }
 1244   switch(typ(x)) {
 1245     case t_INT:
 1246       if (lgefint(K) != 3) return 0;
 1247       return Z_ispowerall(x, itou(K), pt);
 1248     case t_FRAC:
 1249     {
 1250       GEN a = gel(x,1), b = gel(x,2);
 1251       ulong k;
 1252       if (lgefint(K) != 3) return 0;
 1253       k = itou(K);
 1254       if (pt) {
 1255         z = cgetg(3, t_FRAC);
 1256         if (Z_ispowerall(a, k, &a) && Z_ispowerall(b, k, &b)) {
 1257           *pt = z; gel(z,1) = a; gel(z,2) = b; return 1;
 1258         }
 1259         set_avma((pari_sp)(z + 3)); return 0;
 1260       }
 1261       return Z_ispower(a, k) && Z_ispower(b, k);
 1262     }
 1263     case t_INTMOD:
 1264       return Zn_ispower(gel(x,2), gel(x,1), K, pt);
 1265     case t_FFELT:
 1266       return FF_ispower(x, K, pt);
 1267 
 1268     case t_PADIC:
 1269       return Qp_ispower(x, K, pt);
 1270     case t_POLMOD:
 1271       return polmodispower(x, K, pt);
 1272     case t_POL:
 1273       return polispower(x, K, pt);
 1274     case t_RFRAC: {
 1275       GEN a = gel(x,1), b = gel(x,2);
 1276       if (pt) {
 1277         z = cgetg(3, t_RFRAC);
 1278         if (ispower(a, K, &a) && polispower(b, K, &b)) {
 1279           *pt = z; gel(z,1) = a; gel(z,2) = b; return 1;
 1280         }
 1281         set_avma((pari_sp)(z + 3)); return 0;
 1282       }
 1283       return (ispower(a, K, NULL) && polispower(b, K, NULL));
 1284     }
 1285     case t_REAL:
 1286       if (signe(x) < 0 && !mpodd(K)) return 0;
 1287     case t_COMPLEX:
 1288       if (pt) *pt = gsqrtn(x, K, NULL, DEFAULTPREC);
 1289       return 1;
 1290 
 1291     case t_SER:
 1292       if (signe(x) && (!dvdsi(valp(x), K) || !ispower(gel(x,2), K, NULL)))
 1293         return 0;
 1294       if (pt) *pt = gsqrtn(x, K, NULL, DEFAULTPREC);
 1295       return 1;
 1296   }
 1297   pari_err_TYPE("ispower",x);
 1298   return 0; /* LCOV_EXCL_LINE */
 1299 }
 1300 
 1301 long
 1302 gisanypower(GEN x, GEN *pty)
 1303 {
 1304   long tx = typ(x);
 1305   ulong k, h;
 1306   if (tx == t_INT) return Z_isanypower(x, pty);
 1307   if (tx == t_FRAC)
 1308   {
 1309     pari_sp av = avma;
 1310     GEN fa, P, E, a = gel(x,1), b = gel(x,2);
 1311     long i, j, p, e;
 1312     int sw = (abscmpii(a, b) > 0);
 1313 
 1314     if (sw) swap(a, b);
 1315     k = Z_isanypower(a, pty? &a: NULL);
 1316     if (!k)
 1317     { /* a = -1,1 or not a pure power */
 1318       if (!is_pm1(a)) return gc_long(av,0);
 1319       if (signe(a) < 0) b = negi(b);
 1320       k = Z_isanypower(b, pty? &b: NULL);
 1321       if (!k || !pty) return gc_long(av,k);
 1322       *pty = gerepileupto(av, ginv(b));
 1323       return k;
 1324     }
 1325     fa = factoru(k);
 1326     P = gel(fa,1);
 1327     E = gel(fa,2); h = k;
 1328     for (i = lg(P) - 1; i > 0; i--)
 1329     {
 1330       p = P[i];
 1331       e = E[i];
 1332       for (j = 0; j < e; j++)
 1333         if (!is_kth_power(b, p, &b)) break;
 1334       if (j < e) k /= upowuu(p, e - j);
 1335     }
 1336     if (k == 1) return gc_long(av,0);
 1337     if (!pty) return gc_long(av,k);
 1338     if (k != h) a = powiu(a, h/k);
 1339     *pty = gerepilecopy(av, mkfrac(a, b));
 1340     return k;
 1341   }
 1342   pari_err_TYPE("gisanypower", x);
 1343   return 0; /* LCOV_EXCL_LINE */
 1344 }
 1345 
 1346 /* v_p(x) = e != 0 for some p; return ispower(x,,&x), updating x.
 1347  * No need to optimize for 2,3,5,7 powers (done before) */
 1348 static long
 1349 split_exponent(ulong e, GEN *x)
 1350 {
 1351   GEN fa, P, E;
 1352   long i, j, l, k = 1;
 1353   if (e == 1) return 1;
 1354   fa = factoru(e);
 1355   P = gel(fa,1);
 1356   E = gel(fa,2); l = lg(P);
 1357   for (i = 1; i < l; i++)
 1358   {
 1359     ulong p = P[i];
 1360     for (j = 0; j < E[i]; j++)
 1361     {
 1362       GEN y;
 1363       if (!is_kth_power(*x, p, &y)) break;
 1364       k *= p; *x = y;
 1365     }
 1366   }
 1367   return k;
 1368 }
 1369 
 1370 static long
 1371 Z_isanypower_nosmalldiv(GEN *px)
 1372 { /* any prime divisor of x is > 102 */
 1373   const double LOG2_103 = 6.6865; /* lower bound for log_2(103) */
 1374   const double LOG103 = 4.6347; /* lower bound for log(103) */
 1375   forprime_t T;
 1376   ulong mask = 7, e2;
 1377   long k, ex;
 1378   GEN y, x = *px;
 1379 
 1380   k = 1;
 1381   while (Z_issquareall(x, &y)) { k <<= 1; x = y; }
 1382   while ( (ex = is_357_power(x, &y, &mask)) ) { k *= ex; x = y; }
 1383   e2 = (ulong)((expi(x) + 1) / LOG2_103); /* >= log_103 (x) */
 1384   if (u_forprime_init(&T, 11, e2))
 1385   {
 1386     GEN logx = NULL;
 1387     const ulong Q = 30011; /* prime */
 1388     ulong p, xmodQ;
 1389     double dlogx = 0;
 1390     /* cut off at x^(1/p) ~ 2^30 bits which seems to be about optimum;
 1391      * for large p the modular checks are no longer competitively fast */
 1392     while ( (ex = is_pth_power(x, &y, &T, 30)) )
 1393     {
 1394       k *= ex; x = y;
 1395       e2 = (ulong)((expi(x) + 1) / LOG2_103);
 1396       u_forprime_restrict(&T, e2);
 1397     }
 1398     if (DEBUGLEVEL>4)
 1399       err_printf("Z_isanypower: now k=%ld, x=%ld-bit\n", k, expi(x)+1);
 1400     xmodQ = umodiu(x, Q);
 1401     /* test Q | x, just in case */
 1402     if (!xmodQ) { *px = x; return k * split_exponent(Z_lval(x,Q), px); }
 1403     /* x^(1/p) < 2^31 */
 1404     p = T.p;
 1405     if (p <= e2)
 1406     {
 1407       logx = logr_abs( itor(x, DEFAULTPREC) );
 1408       dlogx = rtodbl(logx);
 1409       e2 = (ulong)(dlogx / LOG103); /* >= log_103(x) */
 1410     }
 1411     while (p && p <= e2)
 1412     { /* is x a p-th power ? By computing y = round(x^(1/p)).
 1413        * Check whether y^p = x, first mod Q, then exactly. */
 1414       pari_sp av = avma;
 1415       long e;
 1416       GEN logy = divru(logx, p), y = grndtoi(mpexp(logy), &e);
 1417       ulong ymodQ = umodiu(y,Q);
 1418       if (e >= -10 || Fl_powu(ymodQ, p % (Q-1), Q) != xmodQ
 1419                    || !equalii(powiu(y, p), x)) set_avma(av);
 1420       else
 1421       {
 1422         k *= p; x = y; xmodQ = ymodQ; logx = logy; dlogx /= p;
 1423         e2 = (ulong)(dlogx / LOG103); /* >= log_103(x) */
 1424         u_forprime_restrict(&T, e2);
 1425         continue; /* if success, retry same p */
 1426       }
 1427       p = u_forprime_next(&T);
 1428     }
 1429   }
 1430   *px = x; return k;
 1431 }
 1432 
 1433 static ulong tinyprimes[] = {
 1434   2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71,
 1435   73, 79, 83, 89, 97, 101, 103, 107, 109, 113, 127, 131, 137, 139, 149, 151,
 1436   157, 163, 167, 173, 179, 181, 191, 193, 197, 199
 1437 };
 1438 
 1439 /* disregard the sign of x, caller will take care of x < 0 */
 1440 static long
 1441 Z_isanypower_aux(GEN x, GEN *pty)
 1442 {
 1443   long ex, v, i, l, k;
 1444   GEN y, P, E;
 1445   ulong mask, e = 0;
 1446 
 1447   if (abscmpii(x, gen_2) < 0) return 0; /* -1,0,1 */
 1448 
 1449   if (signe(x) < 0) x = negi(x);
 1450   k = l = 1;
 1451   P = cgetg(26 + 1, t_VECSMALL);
 1452   E = cgetg(26 + 1, t_VECSMALL);
 1453   /* trial division */
 1454   for(i = 0; i < 26; i++)
 1455   {
 1456     ulong p = tinyprimes[i];
 1457     int stop;
 1458     v = Z_lvalrem_stop(&x, p, &stop);
 1459     if (v)
 1460     {
 1461       P[l] = p;
 1462       E[l] = v; l++;
 1463       e = ugcd(e, v); if (e == 1) goto END;
 1464     }
 1465     if (stop) {
 1466       if (is_pm1(x)) k = e;
 1467       goto END;
 1468     }
 1469   }
 1470 
 1471   if (e)
 1472   { /* Bingo. Result divides e */
 1473     long v3, v5, v7;
 1474     ulong e2 = e;
 1475     v = u_lvalrem(e2, 2, &e2);
 1476     if (v)
 1477     {
 1478       for (i = 0; i < v; i++)
 1479       {
 1480         if (!Z_issquareall(x, &y)) break;
 1481         k <<= 1; x = y;
 1482       }
 1483     }
 1484     mask = 0;
 1485     v3 = u_lvalrem(e2, 3, &e2); if (v3) mask = 1;
 1486     v5 = u_lvalrem(e2, 5, &e2); if (v5) mask |= 2;
 1487     v7 = u_lvalrem(e2, 7, &e2); if (v7) mask |= 4;
 1488     while ( (ex = is_357_power(x, &y, &mask)) ) {
 1489       x = y;
 1490       switch(ex)
 1491       {
 1492         case 3: k *= 3; if (--v3 == 0) mask &= ~1; break;
 1493         case 5: k *= 5; if (--v5 == 0) mask &= ~2; break;
 1494         case 7: k *= 7; if (--v7 == 0) mask &= ~4; break;
 1495       }
 1496     }
 1497     k *= split_exponent(e2, &x);
 1498   }
 1499   else
 1500     k = Z_isanypower_nosmalldiv(&x);
 1501 END:
 1502   if (pty && k != 1)
 1503   {
 1504     if (e)
 1505     { /* add missing small factors */
 1506       y = powuu(P[1], E[1] / k);
 1507       for (i = 2; i < l; i++) y = mulii(y, powuu(P[i], E[i] / k));
 1508       x = equali1(x)? y: mulii(x,y);
 1509     }
 1510     *pty = x;
 1511   }
 1512   return k == 1? 0: k;
 1513 }
 1514 
 1515 long
 1516 Z_isanypower(GEN x, GEN *pty)
 1517 {
 1518   pari_sp av = avma;
 1519   long k = Z_isanypower_aux(x, pty);
 1520   if (!k) return gc_long(av,0);
 1521   if (signe(x) < 0)
 1522   {
 1523     long v = vals(k);
 1524     if (v)
 1525     {
 1526       GEN y;
 1527       k >>= v;
 1528       if (k == 1) return gc_long(av,0);
 1529       if (!pty) return gc_long(av,k);
 1530       y = *pty;
 1531       y = powiu(y, 1<<v);
 1532       togglesign(y);
 1533       *pty = gerepileuptoint(av, y);
 1534       return k;
 1535     }
 1536     if (pty) togglesign_safe(pty);
 1537   }
 1538   if (pty) *pty = gerepilecopy(av, *pty); else set_avma(av);
 1539   return k;
 1540 }
 1541 
 1542 /* Faster than */
 1543 /*   !cmpii(n, int2n(vali(n))) */
 1544 /*   !cmpis(shifti(n, -vali(n)), 1) */
 1545 /*   expi(n) == vali(n) */
 1546 /*   hamming(n) == 1 */
 1547 /* even for single-word values, and much faster for multiword values. */
 1548 /* If all you have is a word, you can just use n & !(n & (n-1)). */
 1549 long
 1550 Z_ispow2(GEN n)
 1551 {
 1552   GEN xp;
 1553   long i, lx;
 1554   ulong u;
 1555   if (signe(n) != 1) return 0;
 1556   xp = int_LSW(n);
 1557   lx = lgefint(n);
 1558   u = *xp;
 1559   for (i = 3; i < lx; ++i)
 1560   {
 1561     if (u) return 0;
 1562     xp = int_nextW(xp);
 1563     u = *xp;
 1564   }
 1565   return !(u & (u-1));
 1566 }
 1567 
 1568 static long
 1569 isprimepower_i(GEN n, GEN *pt, long flag)
 1570 {
 1571   pari_sp av = avma;
 1572   long i, v;
 1573 
 1574   if (typ(n) != t_INT) pari_err_TYPE("isprimepower", n);
 1575   if (signe(n) <= 0) return 0;
 1576 
 1577   if (lgefint(n) == 3)
 1578   {
 1579     ulong p;
 1580     v = uisprimepower(n[2], &p);
 1581     if (v)
 1582     {
 1583       if (pt) *pt = utoipos(p);
 1584       return v;
 1585     }
 1586     return 0;
 1587   }
 1588   for (i = 0; i < 26; i++)
 1589   {
 1590     ulong p = tinyprimes[i];
 1591     v = Z_lvalrem(n, p, &n);
 1592     if (v)
 1593     {
 1594       set_avma(av);
 1595       if (!is_pm1(n)) return 0;
 1596       if (pt) *pt = utoipos(p);
 1597       return v;
 1598     }
 1599   }
 1600   /* p | n => p >= 103 */
 1601   v = Z_isanypower_nosmalldiv(&n); /* expensive */
 1602   if (!(flag? isprime(n): BPSW_psp(n))) return gc_long(av,0);
 1603   if (pt) *pt = gerepilecopy(av, n); else set_avma(av);
 1604   return v;
 1605 }
 1606 long
 1607 isprimepower(GEN n, GEN *pt) { return isprimepower_i(n,pt,1); }
 1608 long
 1609 ispseudoprimepower(GEN n, GEN *pt) { return isprimepower_i(n,pt,0); }
 1610 
 1611 long
 1612 uisprimepower(ulong n, ulong *pp)
 1613 { /* We must have CUTOFF^11 >= ULONG_MAX and CUTOFF^3 < ULONG_MAX.
 1614    * Tests suggest that 200-300 is the best range for 64-bit platforms. */
 1615   const ulong CUTOFF = 200UL;
 1616   const long TINYCUTOFF = 46;  /* tinyprimes[45] = 199 */
 1617   const ulong CUTOFF3 = CUTOFF*CUTOFF*CUTOFF;
 1618 #ifdef LONG_IS_64BIT
 1619   /* primes preceeding the appropriate root of ULONG_MAX. */
 1620   const ulong ROOT9 = 137;
 1621   const ulong ROOT8 = 251;
 1622   const ulong ROOT7 = 563;
 1623   const ulong ROOT5 = 7129;
 1624   const ulong ROOT4 = 65521;
 1625 #else
 1626   const ulong ROOT9 = 11;
 1627   const ulong ROOT8 = 13;
 1628   const ulong ROOT7 = 23;
 1629   const ulong ROOT5 = 83;
 1630   const ulong ROOT4 = 251;
 1631 #endif
 1632   ulong mask;
 1633   long v, i;
 1634   int e;
 1635   if (n < 2) return 0;
 1636   if (!odd(n)) {
 1637     if (n & (n-1)) return 0;
 1638     *pp = 2; return vals(n);
 1639   }
 1640   if (n < 8) { *pp = n; return 1; } /* 3,5,7 */
 1641   for (i = 1/*skip p=2*/; i < TINYCUTOFF; i++)
 1642   {
 1643     ulong p = tinyprimes[i];
 1644     if (n % p == 0)
 1645     {
 1646       v = u_lvalrem(n, p, &n);
 1647       if (n == 1) { *pp = p; return v; }
 1648       return 0;
 1649     }
 1650   }
 1651   /* p | n => p >= CUTOFF */
 1652 
 1653   if (n < CUTOFF3)
 1654   {
 1655     if (n < CUTOFF*CUTOFF || uisprime_101(n)) { *pp = n; return 1; }
 1656     if (uissquareall(n, &n)) { *pp = n; return 2; }
 1657     return 0;
 1658   }
 1659 
 1660   /* Check for squares, fourth powers, and eighth powers as appropriate. */
 1661   v = 1;
 1662   if (uissquareall(n, &n)) {
 1663     v <<= 1;
 1664     if (CUTOFF <= ROOT4 && uissquareall(n, &n)) {
 1665       v <<= 1;
 1666       if (CUTOFF <= ROOT8 && uissquareall(n, &n)) v <<= 1;
 1667     }
 1668   }
 1669 
 1670   if (CUTOFF > ROOT5) mask = 1;
 1671   else
 1672   {
 1673     const ulong CUTOFF5 = CUTOFF3*CUTOFF*CUTOFF;
 1674     if (n < CUTOFF5) mask = 1; else mask = 3;
 1675     if (CUTOFF <= ROOT7)
 1676     {
 1677       const ulong CUTOFF7 = CUTOFF5*CUTOFF*CUTOFF;
 1678       if (n >= CUTOFF7) mask = 7;
 1679     }
 1680   }
 1681 
 1682   if (CUTOFF <= ROOT9 && (e = uis_357_power(n, &n, &mask))) { v *= e; mask=1; }
 1683   if ((e = uis_357_power(n, &n, &mask))) v *= e;
 1684 
 1685   if (uisprime_101(n)) { *pp = n; return v; }
 1686   return 0;
 1687 }
 1688 
 1689 /*********************************************************************/
 1690 /**                                                                 **/
 1691 /**                        KRONECKER SYMBOL                         **/
 1692 /**                                                                 **/
 1693 /*********************************************************************/
 1694 /* t = 3,5 mod 8 ?  (= 2 not a square mod t) */
 1695 static int
 1696 ome(long t)
 1697 {
 1698   switch(t & 7)
 1699   {
 1700     case 3:
 1701     case 5: return 1;
 1702     default: return 0;
 1703   }
 1704 }
 1705 /* t a t_INT, is t = 3,5 mod 8 ? */
 1706 static int
 1707 gome(GEN t)
 1708 { return signe(t)? ome( mod2BIL(t) ): 0; }
 1709 
 1710 /* assume y odd, return kronecker(x,y) * s */
 1711 static long
 1712 krouu_s(ulong x, ulong y, long s)
 1713 {
 1714   ulong x1 = x, y1 = y, z;
 1715   while (x1)
 1716   {
 1717     long r = vals(x1);
 1718     if (r)
 1719     {
 1720       if (odd(r) && ome(y1)) s = -s;
 1721       x1 >>= r;
 1722     }
 1723     if (x1 & y1 & 2) s = -s;
 1724     z = y1 % x1; y1 = x1; x1 = z;
 1725   }
 1726   return (y1 == 1)? s: 0;
 1727 }
 1728 
 1729 long
 1730 kronecker(GEN x, GEN y)
 1731 {
 1732   pari_sp av = avma;
 1733   long s = 1, r;
 1734   ulong xu;
 1735 
 1736   if (typ(x) != t_INT) pari_err_TYPE("kronecker",x);
 1737   if (typ(y) != t_INT) pari_err_TYPE("kronecker",y);
 1738   switch (signe(y))
 1739   {
 1740     case -1: y = negi(y); if (signe(x) < 0) s = -1; break;
 1741     case 0: return is_pm1(x);
 1742   }
 1743   r = vali(y);
 1744   if (r)
 1745   {
 1746     if (!mpodd(x)) return gc_long(av,0);
 1747     if (odd(r) && gome(x)) s = -s;
 1748     y = shifti(y,-r);
 1749   }
 1750   x = modii(x,y);
 1751   while (lgefint(x) > 3) /* x < y */
 1752   {
 1753     GEN z;
 1754     r = vali(x);
 1755     if (r)
 1756     {
 1757       if (odd(r) && gome(y)) s = -s;
 1758       x = shifti(x,-r);
 1759     }
 1760     /* x=3 mod 4 && y=3 mod 4 ? (both are odd here) */
 1761     if (mod2BIL(x) & mod2BIL(y) & 2) s = -s;
 1762     z = remii(y,x); y = x; x = z;
 1763     if (gc_needed(av,2))
 1764     {
 1765       if(DEBUGMEM>1) pari_warn(warnmem,"kronecker");
 1766       gerepileall(av, 2, &x, &y);
 1767     }
 1768   }
 1769   xu = itou(x);
 1770   if (!xu) return is_pm1(y)? s: 0;
 1771   r = vals(xu);
 1772   if (r)
 1773   {
 1774     if (odd(r) && gome(y)) s = -s;
 1775     xu >>= r;
 1776   }
 1777   /* x=3 mod 4 && y=3 mod 4 ? (both are odd here) */
 1778   if (xu & mod2BIL(y) & 2) s = -s;
 1779   return gc_long(av, krouu_s(umodiu(y,xu), xu, s));
 1780 }
 1781 
 1782 long
 1783 krois(GEN x, long y)
 1784 {
 1785   ulong yu;
 1786   long s = 1;
 1787 
 1788   if (y <= 0)
 1789   {
 1790     if (y == 0) return is_pm1(x);
 1791     yu = (ulong)-y; if (signe(x) < 0) s = -1;
 1792   }
 1793   else
 1794     yu = (ulong)y;
 1795   if (!odd(yu))
 1796   {
 1797     long r;
 1798     if (!mpodd(x)) return 0;
 1799     r = vals(yu); yu >>= r;
 1800     if (odd(r) && gome(x)) s = -s;
 1801   }
 1802   return krouu_s(umodiu(x, yu), yu, s);
 1803 }
 1804 /* assume y != 0 */
 1805 long
 1806 kroiu(GEN x, ulong y)
 1807 {
 1808   long r;
 1809   if (odd(y)) return krouu_s(umodiu(x,y), y, 1);
 1810   if (!mpodd(x)) return 0;
 1811   r = vals(y); y >>= r;
 1812   return krouu_s(umodiu(x,y), y, (odd(r) && gome(x))? -1: 1);
 1813 }
 1814 
 1815 /* assume y > 0, odd, return s * kronecker(x,y) */
 1816 static long
 1817 krouodd(ulong x, GEN y, long s)
 1818 {
 1819   long r;
 1820   if (lgefint(y) == 3) return krouu_s(x, y[2], s);
 1821   if (!x) return 0; /* y != 1 */
 1822   r = vals(x);
 1823   if (r)
 1824   {
 1825     if (odd(r) && gome(y)) s = -s;
 1826     x >>= r;
 1827   }
 1828   /* x=3 mod 4 && y=3 mod 4 ? (both are odd here) */
 1829   if (x & mod2BIL(y) & 2) s = -s;
 1830   return krouu_s(umodiu(y,x), x, s);
 1831 }
 1832 
 1833 long
 1834 krosi(long x, GEN y)
 1835 {
 1836   const pari_sp av = avma;
 1837   long s = 1, r;
 1838   switch (signe(y))
 1839   {
 1840     case -1: y = negi(y); if (x < 0) s = -1; break;
 1841     case 0: return (x==1 || x==-1);
 1842   }
 1843   r = vali(y);
 1844   if (r)
 1845   {
 1846     if (!odd(x)) return gc_long(av,0);
 1847     if (odd(r) && ome(x)) s = -s;
 1848     y = shifti(y,-r);
 1849   }
 1850   if (x < 0) { x = -x; if (mod4(y) == 3) s = -s; }
 1851   return gc_long(av, krouodd((ulong)x, y, s));
 1852 }
 1853 
 1854 long
 1855 kroui(ulong x, GEN y)
 1856 {
 1857   const pari_sp av = avma;
 1858   long s = 1, r;
 1859   switch (signe(y))
 1860   {
 1861     case -1: y = negi(y); break;
 1862     case 0: return x==1UL;
 1863   }
 1864   r = vali(y);
 1865   if (r)
 1866   {
 1867     if (!odd(x)) return gc_long(av,0);
 1868     if (odd(r) && ome(x)) s = -s;
 1869     y = shifti(y,-r);
 1870   }
 1871   return gc_long(av,  krouodd(x, y, s));
 1872 }
 1873 
 1874 long
 1875 kross(long x, long y)
 1876 {
 1877   ulong yu;
 1878   long s = 1;
 1879 
 1880   if (y <= 0)
 1881   {
 1882     if (y == 0) return (labs(x)==1);
 1883     yu = (ulong)-y; if (x < 0) s = -1;
 1884   }
 1885   else
 1886     yu = (ulong)y;
 1887   if (!odd(yu))
 1888   {
 1889     long r;
 1890     if (!odd(x)) return 0;
 1891     r = vals(yu); yu >>= r;
 1892     if (odd(r) && ome(x)) s = -s;
 1893   }
 1894   x %= (long)yu; if (x < 0) x += yu;
 1895   return krouu_s((ulong)x, yu, s);
 1896 }
 1897 
 1898 long
 1899 krouu(ulong x, ulong y)
 1900 {
 1901   long r;
 1902   if (odd(y)) return krouu_s(x, y, 1);
 1903   if (!odd(x)) return 0;
 1904   r = vals(y); y >>= r;
 1905   return krouu_s(x, y, (odd(r) && ome(x))? -1: 1);
 1906 }
 1907 
 1908 /*********************************************************************/
 1909 /**                                                                 **/
 1910 /**                          HILBERT SYMBOL                         **/
 1911 /**                                                                 **/
 1912 /*********************************************************************/
 1913 /* x,y are t_INT or t_REAL */
 1914 static long
 1915 mphilbertoo(GEN x, GEN y)
 1916 {
 1917   long sx = signe(x), sy = signe(y);
 1918   if (!sx || !sy) return 0;
 1919   return (sx < 0 && sy < 0)? -1: 1;
 1920 }
 1921 
 1922 long
 1923 hilbertii(GEN x, GEN y, GEN p)
 1924 {
 1925   pari_sp av;
 1926   long oddvx, oddvy, z;
 1927 
 1928   if (!p) return mphilbertoo(x,y);
 1929   if (is_pm1(p) || signe(p) < 0) pari_err_PRIME("hilbertii",p);
 1930   if (!signe(x) || !signe(y)) return 0;
 1931   av = avma;
 1932   oddvx = odd(Z_pvalrem(x,p,&x));
 1933   oddvy = odd(Z_pvalrem(y,p,&y));
 1934   /* x, y are p-units, compute hilbert(x * p^oddvx, y * p^oddvy, p) */
 1935   if (absequaliu(p, 2))
 1936   {
 1937     z = (Mod4(x) == 3 && Mod4(y) == 3)? -1: 1;
 1938     if (oddvx && gome(y)) z = -z;
 1939     if (oddvy && gome(x)) z = -z;
 1940   }
 1941   else
 1942   {
 1943     z = (oddvx && oddvy && mod4(p) == 3)? -1: 1;
 1944     if (oddvx && kronecker(y,p) < 0) z = -z;
 1945     if (oddvy && kronecker(x,p) < 0) z = -z;
 1946   }
 1947   return gc_long(av, z);
 1948 }
 1949 
 1950 static void
 1951 err_prec(void) { pari_err_PREC("hilbert"); }
 1952 static void
 1953 err_p(GEN p, GEN q) { pari_err_MODULUS("hilbert", p,q); }
 1954 static void
 1955 err_oo(GEN p) { pari_err_MODULUS("hilbert", p, strtoGENstr("oo")); }
 1956 
 1957 /* x t_INTMOD, *pp = prime or NULL [ unset, set it to x.mod ].
 1958  * Return lift(x) provided it's p-adic accuracy is large enough to decide
 1959  * hilbert()'s value [ problem at p = 2 ] */
 1960 static GEN
 1961 lift_intmod(GEN x, GEN *pp)
 1962 {
 1963   GEN p = *pp, N = gel(x,1);
 1964   x = gel(x,2);
 1965   if (!p)
 1966   {
 1967     *pp = p = N;
 1968     switch(itos_or_0(p))
 1969     {
 1970       case 2:
 1971       case 4: err_prec();
 1972     }
 1973     return x;
 1974   }
 1975   if (!signe(p)) err_oo(N);
 1976   if (absequaliu(p,2))
 1977   { if (vali(N) <= 2) err_prec(); }
 1978   else
 1979   { if (!dvdii(N,p)) err_p(N,p); }
 1980   if (!signe(x)) err_prec();
 1981   return x;
 1982 }
 1983 /* x t_PADIC, *pp = prime or NULL [ unset, set it to x.p ].
 1984  * Return lift(x)*p^(v(x) mod 2) provided it's p-adic accuracy is large enough
 1985  * to decide hilbert()'s value [ problem at p = 2 ]*/
 1986 static GEN
 1987 lift_padic(GEN x, GEN *pp)
 1988 {
 1989   GEN p = *pp, q = gel(x,2), y = gel(x,4);
 1990   if (!p) *pp = p = q;
 1991   else if (!equalii(p,q)) err_p(p, q);
 1992   if (absequaliu(p,2) && precp(x) <= 2) err_prec();
 1993   if (!signe(y)) err_prec();
 1994   return odd(valp(x))? mulii(p,y): y;
 1995 }
 1996 
 1997 long
 1998 hilbert(GEN x, GEN y, GEN p)
 1999 {
 2000   pari_sp av = avma;
 2001   long tx = typ(x), ty = typ(y);
 2002 
 2003   if (p && typ(p) != t_INT) pari_err_TYPE("hilbert",p);
 2004   if (tx == t_REAL)
 2005   {
 2006     if (p && signe(p)) err_oo(p);
 2007     switch (ty)
 2008     {
 2009       case t_INT:
 2010       case t_REAL: return mphilbertoo(x,y);
 2011       case t_FRAC: return mphilbertoo(x,gel(y,1));
 2012       default: pari_err_TYPE2("hilbert",x,y);
 2013     }
 2014   }
 2015   if (ty == t_REAL)
 2016   {
 2017     if (p && signe(p)) err_oo(p);
 2018     switch (tx)
 2019     {
 2020       case t_INT:
 2021       case t_REAL: return mphilbertoo(x,y);
 2022       case t_FRAC: return mphilbertoo(gel(x,1),y);
 2023       default: pari_err_TYPE2("hilbert",x,y);
 2024     }
 2025   }
 2026   if (tx == t_INTMOD) { x = lift_intmod(x, &p); tx = t_INT; }
 2027   if (ty == t_INTMOD) { y = lift_intmod(y, &p); ty = t_INT; }
 2028 
 2029   if (tx == t_PADIC) { x = lift_padic(x, &p); tx = t_INT; }
 2030   if (ty == t_PADIC) { y = lift_padic(y, &p); ty = t_INT; }
 2031 
 2032   if (tx == t_FRAC) { tx = t_INT; x = p? mulii(gel(x,1),gel(x,2)): gel(x,1); }
 2033   if (ty == t_FRAC) { ty = t_INT; y = p? mulii(gel(y,1),gel(y,2)): gel(y,1); }
 2034 
 2035   if (tx != t_INT || ty != t_INT) pari_err_TYPE2("hilbert",x,y);
 2036   if (p && !signe(p)) p = NULL;
 2037   return gc_long(av, hilbertii(x,y,p));
 2038 }
 2039 
 2040 /*******************************************************************/
 2041 /*                                                                 */
 2042 /*                       SQUARE ROOT MODULO p                      */
 2043 /*                                                                 */
 2044 /*******************************************************************/
 2045 static void
 2046 checkp(ulong q, ulong p)
 2047 { if (!q) pari_err_PRIME("Fl_nonsquare",utoipos(p)); }
 2048 /* p = 1 (mod 4) prime, return the first quadratic nonresidue, a prime */
 2049 static ulong
 2050 nonsquare1_Fl(ulong p)
 2051 {
 2052   forprime_t S;
 2053   ulong q;
 2054   if ((p & 7UL) != 1) return 2UL;
 2055   q = p % 3; if (q == 2) return 3UL;
 2056   checkp(q, p);
 2057   q = p % 5; if (q == 2 || q == 3) return 5UL;
 2058   checkp(q, p);
 2059   q = p % 7; if (q != 4 && q >= 3) return 7UL;
 2060   checkp(q, p);
 2061   u_forprime_init(&S, 11, p);
 2062   while ((q = u_forprime_next(&S)))
 2063   {
 2064     long i = krouu(q, p);
 2065     if (i < 0) return q;
 2066     checkp(q, p);
 2067   }
 2068   checkp(0, p);
 2069   return 0; /*LCOV_EXCL_LINE*/
 2070 }
 2071 /* p > 2 a prime */
 2072 ulong
 2073 nonsquare_Fl(ulong p)
 2074 { return ((p & 3UL) == 3)? p-1: nonsquare1_Fl(p); }
 2075 
 2076 ulong
 2077 Fl_2gener_pre(ulong p, ulong pi)
 2078 {
 2079   ulong p1 = p-1;
 2080   long e = vals(p1);
 2081   if (e == 1) return p1;
 2082   return Fl_powu_pre(nonsquare1_Fl(p), p1 >> e, p, pi);
 2083 }
 2084 
 2085 /* Tonelli-Shanks. Assume p is prime and (a,p) != -1. */
 2086 ulong
 2087 Fl_sqrt_pre_i(ulong a, ulong y, ulong p, ulong pi)
 2088 {
 2089   long i, e, k;
 2090   ulong p1, q, v, w;
 2091 
 2092   if (!a) return 0;
 2093   p1 = p - 1; e = vals(p1);
 2094   if (e == 0) /* p = 2 */
 2095   {
 2096     if (p != 2) pari_err_PRIME("Fl_sqrt [modulus]",utoi(p));
 2097     return ((a & 1) == 0)? 0: 1;
 2098   }
 2099   if (e == 1)
 2100   {
 2101     v = Fl_powu_pre(a, (p+1) >> 2, p, pi);
 2102     if (Fl_sqr_pre(v, p, pi) != a) return ~0UL;
 2103     p1 = p - v; if (v > p1) v = p1;
 2104     return v;
 2105   }
 2106   q = p1 >> e; /* q = (p-1)/2^oo is odd */
 2107   p1 = Fl_powu_pre(a, q >> 1, p, pi); /* a ^ [(q-1)/2] */
 2108   if (!p1) return 0;
 2109   v = Fl_mul_pre(a, p1, p, pi);
 2110   w = Fl_mul_pre(v, p1, p, pi);
 2111   if (!y) y = Fl_powu_pre(nonsquare1_Fl(p), q, p, pi);
 2112   while (w != 1)
 2113   { /* a*w = v^2, y primitive 2^e-th root of 1
 2114        a square --> w even power of y, hence w^(2^(e-1)) = 1 */
 2115     p1 = Fl_sqr_pre(w,p,pi);
 2116     for (k=1; p1 != 1 && k < e; k++) p1 = Fl_sqr_pre(p1,p,pi);
 2117     if (k == e) return ~0UL;
 2118     /* w ^ (2^k) = 1 --> w = y ^ (u * 2^(e-k)), u odd */
 2119     p1 = y;
 2120     for (i=1; i < e-k; i++) p1 = Fl_sqr_pre(p1, p, pi);
 2121     y = Fl_sqr_pre(p1, p, pi); e = k;
 2122     w = Fl_mul_pre(y, w, p, pi);
 2123     v = Fl_mul_pre(v, p1, p, pi);
 2124   }
 2125   p1 = p - v; if (v > p1) v = p1;
 2126   return v;
 2127 }
 2128 
 2129 ulong
 2130 Fl_sqrt(ulong a, ulong p)
 2131 {
 2132   ulong pi = get_Fl_red(p);
 2133   return Fl_sqrt_pre_i(a, 0, p, pi);
 2134 }
 2135 
 2136 ulong
 2137 Fl_sqrt_pre(ulong a, ulong p, ulong pi)
 2138 {
 2139   return Fl_sqrt_pre_i(a, 0, p, pi);
 2140 }
 2141 
 2142 static ulong
 2143 Fl_lgener_pre_all(ulong l, long e, ulong r, ulong p, ulong pi, ulong *pt_m)
 2144 {
 2145   ulong x, y, m;
 2146   ulong le1 = upowuu(l, e-1);
 2147   for (x = 2; ; x++)
 2148   {
 2149     y = Fl_powu_pre(x, r, p, pi);
 2150     if (y==1) continue;
 2151     m = Fl_powu_pre(y, le1, p, pi);
 2152     if (m != 1) break;
 2153   }
 2154   *pt_m = m;
 2155   return y;
 2156 }
 2157 
 2158 /* solve x^l = a , l prime in G of order q.
 2159  *
 2160  * q =  (l^e)*r, e >= 1, (r,l) = 1
 2161  * y generates the l-Sylow of G
 2162  * m = y^(l^(e-1)) != 1 */
 2163 static ulong
 2164 Fl_sqrtl_raw(ulong a, ulong l, ulong e, ulong r, ulong p, ulong pi, ulong y, ulong m)
 2165 {
 2166   ulong p1, v, w, z, dl;
 2167   ulong u2;
 2168   if (a==0) return a;
 2169   u2 = Fl_inv(l%r, r);
 2170   v = Fl_powu_pre(a, u2, p, pi);
 2171   w = Fl_powu_pre(v, l, p, pi);
 2172   w = Fl_mul_pre(w, Fl_inv(a, p), p, pi);
 2173   if (w==1) return v;
 2174   if (y==0) y = Fl_lgener_pre_all(l, e, r, p, pi, &m);
 2175   while (w!=1)
 2176   {
 2177     ulong k = 0;
 2178     p1 = w;
 2179     do
 2180     {
 2181       z = p1; p1 = Fl_powu_pre(p1, l, p, pi);
 2182       if (++k == e) return ULONG_MAX;
 2183     } while (p1!=1);
 2184     dl = Fl_log_pre(z, m, l, p, pi);
 2185     dl = Fl_neg(dl, l);
 2186     p1 = Fl_powu_pre(y,dl*upowuu(l,e-k-1),p,pi);
 2187     m = Fl_powu_pre(m, dl, p, pi);
 2188     e = k;
 2189     v = Fl_mul_pre(p1,v,p,pi);
 2190     y = Fl_powu_pre(p1,l,p,pi);
 2191     w = Fl_mul_pre(y,w,p,pi);
 2192   }
 2193   return v;
 2194 }
 2195 
 2196 static ulong
 2197 Fl_sqrtl_i(ulong a, ulong l, ulong p, ulong pi, ulong y, ulong m)
 2198 {
 2199   ulong r, e = u_lvalrem(p-1, l, &r);
 2200   return Fl_sqrtl_raw(a, l, e, r, p, pi, y, m);
 2201 }
 2202 
 2203 ulong
 2204 Fl_sqrtl_pre(ulong a, ulong l, ulong p, ulong pi)
 2205 {
 2206   return Fl_sqrtl_i(a, l, p, pi, 0, 0);
 2207 }
 2208 
 2209 ulong
 2210 Fl_sqrtl(ulong a, ulong l, ulong p)
 2211 {
 2212   ulong pi = get_Fl_red(p);
 2213   return Fl_sqrtl_i(a, l, p, pi, 0, 0);
 2214 }
 2215 
 2216 ulong
 2217 Fl_sqrtn_pre(ulong a, long n, ulong p, ulong pi, ulong *zetan)
 2218 {
 2219   ulong m, q = p-1, z;
 2220   ulong nn = n >= 0 ? (ulong)n: -(ulong)n;
 2221   if (a==0)
 2222   {
 2223     if (n < 0) pari_err_INV("Fl_sqrtn", mkintmod(gen_0,utoi(p)));
 2224     if (zetan) *zetan = 1UL;
 2225     return 0;
 2226   }
 2227   if (n==1)
 2228   {
 2229     if (zetan) *zetan = 1;
 2230     return n < 0? Fl_inv(a,p): a;
 2231   }
 2232   if (n==2)
 2233   {
 2234     if (zetan) *zetan = p-1;
 2235     return Fl_sqrt_pre_i(a, 0, p, pi);
 2236   }
 2237   if (a == 1 && !zetan) return a;
 2238   m = ugcd(nn, q);
 2239   z = 1;
 2240   if (m!=1)
 2241   {
 2242     GEN F = factoru(m);
 2243     long i, j, e;
 2244     ulong r, zeta, y, l;
 2245     for (i = nbrows(F); i; i--)
 2246     {
 2247       l = ucoeff(F,i,1);
 2248       j = ucoeff(F,i,2);
 2249       e = u_lvalrem(q,l, &r);
 2250       y = Fl_lgener_pre_all(l, e, r, p, pi, &zeta);
 2251       if (zetan)
 2252         z = Fl_mul_pre(z, Fl_powu_pre(y, upowuu(l,e-j), p, pi), p, pi);
 2253       if (a!=1)
 2254         do
 2255         {
 2256           a = Fl_sqrtl_raw(a, l, e, r, p, pi, y, zeta);
 2257           if (a==ULONG_MAX) return ULONG_MAX;
 2258         } while (--j);
 2259     }
 2260   }
 2261   if (m != nn)
 2262   {
 2263     ulong qm = q/m, nm = nn/m;
 2264     a = Fl_powu_pre(a, Fl_inv(nm%qm, qm), p, pi);
 2265   }
 2266   if (n < 0) a = Fl_inv(a, p);
 2267   if (zetan) *zetan = z;
 2268   return a;
 2269 }
 2270 
 2271 ulong
 2272 Fl_sqrtn(ulong a, long n, ulong p, ulong *zetan)
 2273 {
 2274   ulong pi = get_Fl_red(p);
 2275   return Fl_sqrtn_pre(a, n, p, pi, zetan);
 2276 }
 2277 
 2278 /* Cipolla is better than Tonelli-Shanks when e = v_2(p-1) is "too big".
 2279  * Otherwise, is a constant times worse; for p = 3 (mod 4), is about 3 times worse,
 2280  * and in average is about 2 or 2.5 times worse. But try both algorithms for
 2281  * S(n) = (2^n+3)^2-8 with n = 750, 771, 779, 790, 874, 1176, 1728, 2604, etc.
 2282  *
 2283  * If X^2 := t^2 - a  is not a square in F_p (so X is in F_p^2), then
 2284  *   (t+X)^(p+1) = (t-X)(t+X) = a,   hence  sqrt(a) = (t+X)^((p+1)/2)  in F_p^2.
 2285  * If (a|p)=1, then sqrt(a) is in F_p.
 2286  * cf: LNCS 2286, pp 430-434 (2002)  [Gonzalo Tornaria] */
 2287 
 2288 /* compute y^2, y = y[1] + y[2] X */
 2289 static GEN
 2290 sqrt_Cipolla_sqr(void *data, GEN y)
 2291 {
 2292   GEN u = gel(y,1), v = gel(y,2), p = gel(data,2), n = gel(data,3);
 2293   GEN u2 = sqri(u), v2 = sqri(v);
 2294   v = subii(sqri(addii(v,u)), addii(u2,v2));
 2295   u = addii(u2, mulii(v2,n));
 2296   /* NOT mkvec2: must be gerepileupto-able */
 2297   retmkvec2(modii(u,p), modii(v,p));
 2298 }
 2299 /* compute (t+X) y^2 */
 2300 static GEN
 2301 sqrt_Cipolla_msqr(void *data, GEN y)
 2302 {
 2303   GEN u = gel(y,1), v = gel(y,2), a = gel(data,1), p = gel(data,2), gt = gel(data,4);
 2304   ulong t = gt[2];
 2305   GEN d = addii(u, mului(t,v)), d2= sqri(d);
 2306   GEN b = remii(mulii(a,v), p);
 2307   u = subii(mului(t,d2), mulii(b,addii(u,d)));
 2308   v = subii(d2, mulii(b,v));
 2309   /* NOT mkvec2: must be gerepileupto-able */
 2310   retmkvec2(modii(u,p), modii(v,p));
 2311 }
 2312 /* assume a reduced mod p [ otherwise correct but inefficient ] */
 2313 static GEN
 2314 sqrt_Cipolla(GEN a, GEN p)
 2315 {
 2316   pari_sp av1;
 2317   GEN u, v, n, y, pov2;
 2318   ulong t;
 2319 
 2320   if (kronecker(a, p) < 0) return NULL;
 2321   pov2 = shifti(p,-1);
 2322   if (cmpii(a,pov2) > 0) a = subii(a,p); /* center: avoid multiplying by huge base*/
 2323 
 2324   av1 = avma;
 2325   for(t=1; ; t++)
 2326   {
 2327     n = subsi((long)(t*t), a);
 2328     if (kronecker(n, p) < 0) break;
 2329     set_avma(av1);
 2330   }
 2331 
 2332   /* compute (t+X)^((p-1)/2) =: u+vX */
 2333   u = utoipos(t);
 2334   y = gen_pow_fold(mkvec2(u, gen_1), pov2, mkvec4(a,p,n,u),
 2335                          sqrt_Cipolla_sqr, sqrt_Cipolla_msqr);
 2336   /* Now u+vX = (t+X)^((p-1)/2); thus
 2337    *   (u+vX)(t+X) = sqrt(a) + 0 X
 2338    * Whence,
 2339    *   sqrt(a) = (u+vt)t - v*a
 2340    *   0       = (u+vt)
 2341    * Thus a square root is v*a */
 2342 
 2343   v = Fp_mul(gel(y, 2), a, p);
 2344   if (cmpii(v,pov2) > 0) v = subii(p,v);
 2345   return v;
 2346 }
 2347 
 2348 /* Return NULL if p is found to be composite */
 2349 static GEN
 2350 Fp_2gener_all(long e, GEN p)
 2351 {
 2352   GEN y, m;
 2353   long k;
 2354   GEN q = shifti(subiu(p,1), -e); /* q = (p-1)/2^oo is odd */
 2355   if (e==0 && !equaliu(p,2)) return NULL;
 2356   for (k=2; ; k++)
 2357   {
 2358     long i = krosi(k, p);
 2359     if (i >= 0)
 2360     {
 2361       if (i) continue;
 2362       return NULL;
 2363     }
 2364     y = m = Fp_pow(utoi(k), q, p);
 2365     for (i=1; i<e; i++)
 2366       if (equali1(m = Fp_sqr(m, p))) break;
 2367     if (i == e) break; /* success */
 2368   }
 2369   return y;
 2370 }
 2371 
 2372 /* Return NULL if p is found to be composite */
 2373 GEN
 2374 Fp_2gener(GEN p)
 2375 { return Fp_2gener_all(vali(subis(p,1)),p); }
 2376 
 2377 /* smallest square root */
 2378 static GEN
 2379 choose_sqrt(GEN v, GEN p)
 2380 {
 2381   pari_sp av = avma;
 2382   GEN q = subii(p,v);
 2383   if (cmpii(v,q) > 0) v = q; else set_avma(av);
 2384   return v;
 2385 }
 2386 /* Tonelli-Shanks. Assume p is prime and return NULL if (a,p) = -1. */
 2387 GEN
 2388 Fp_sqrt_i(GEN a, GEN y, GEN p)
 2389 {
 2390   pari_sp av = avma;
 2391   long i, k, e;
 2392   GEN p1, q, v, w;
 2393 
 2394   if (typ(a) != t_INT) pari_err_TYPE("Fp_sqrt",a);
 2395   if (typ(p) != t_INT) pari_err_TYPE("Fp_sqrt",p);
 2396   if (signe(p) <= 0 || equali1(p)) pari_err_PRIME("Fp_sqrt",p);
 2397   if (lgefint(p) == 3)
 2398   {
 2399     ulong pp = uel(p,2), u = Fl_sqrt(umodiu(a, pp), pp);
 2400     if (u == ~0UL) return NULL;
 2401     return utoi(u);
 2402   }
 2403 
 2404   a = modii(a, p); if (!signe(a)) { set_avma(av); return gen_0; }
 2405   p1 = subiu(p,1); e = vali(p1);
 2406   if (e <= 2)
 2407   { /* direct formulas more efficient */
 2408     pari_sp av2;
 2409     if (e == 0) pari_err_PRIME("Fp_sqrt [modulus]",p); /* p != 2 */
 2410     if (e == 1)
 2411     {
 2412       q = addiu(shifti(p1,-2),1); /* (p+1) / 4 */
 2413       v = Fp_pow(a, q, p);
 2414     }
 2415     else
 2416     { /* Atkin's formula */
 2417       GEN i, a2 = shifti(a,1);
 2418       if (cmpii(a2,p) >= 0) a2 = subii(a2,p);
 2419       q = shifti(p1, -3); /* (p-5)/8 */
 2420       v = Fp_pow(a2, q, p);
 2421       i = Fp_mul(a2, Fp_sqr(v,p), p); /* i^2 = -1 */
 2422       v = Fp_mul(a, Fp_mul(v, subiu(i,1), p), p);
 2423     }
 2424     av2 = avma;
 2425     /* must check equality in case (a/p) = -1 or p not prime */
 2426     e = equalii(Fp_sqr(v,p), a); set_avma(av2);
 2427     return e? gerepileuptoint(av,choose_sqrt(v,p)): NULL;
 2428   }
 2429   /* On average, Cipolla is better than Tonelli/Shanks if and only if
 2430    * e(e-1) > 8*log2(n)+20, see LNCS 2286 pp 430 [GTL] */
 2431   if (e*(e-1) > 20 + 8 * expi(p))
 2432   {
 2433     v = sqrt_Cipolla(a,p); if (!v) return gc_NULL(av);
 2434     return gerepileuptoint(av,v);
 2435   }
 2436   if (!y)
 2437   {
 2438     y = Fp_2gener_all(e, p);
 2439     if (!y) pari_err_PRIME("Fp_sqrt [modulus]",p);
 2440   }
 2441   q = shifti(p1,-e); /* q = (p-1)/2^oo is odd */
 2442   p1 = Fp_pow(a, shifti(q,-1), p); /* a ^ (q-1)/2 */
 2443   v = Fp_mul(a, p1, p);
 2444   w = Fp_mul(v, p1, p);
 2445   while (!equali1(w))
 2446   { /* a*w = v^2, y primitive 2^e-th root of 1
 2447        a square --> w even power of y, hence w^(2^(e-1)) = 1 */
 2448     p1 = Fp_sqr(w,p);
 2449     for (k=1; !equali1(p1) && k < e; k++) p1 = Fp_sqr(p1,p);
 2450     if (k == e) return gc_NULL(av); /* p composite or (a/p) != 1 */
 2451     /* w ^ (2^k) = 1 --> w = y ^ (u * 2^(e-k)), u odd */
 2452     p1 = y;
 2453     for (i=1; i < e-k; i++) p1 = Fp_sqr(p1,p);
 2454     y = Fp_sqr(p1, p); e = k;
 2455     w = Fp_mul(y, w, p);
 2456     v = Fp_mul(v, p1, p);
 2457     if (gc_needed(av,1))
 2458     {
 2459       if(DEBUGMEM>1) pari_warn(warnmem,"Fp_sqrt");
 2460       gerepileall(av,3, &y,&w,&v);
 2461     }
 2462   }
 2463   return gerepileuptoint(av, choose_sqrt(v,p));
 2464 }
 2465 
 2466 GEN
 2467 Fp_sqrt(GEN a, GEN p)
 2468 {
 2469   return Fp_sqrt_i(a, NULL, p);
 2470 }
 2471 
 2472 /*********************************************************************/
 2473 /**                                                                 **/
 2474 /**                        GCD & BEZOUT                             **/
 2475 /**                                                                 **/
 2476 /*********************************************************************/
 2477 
 2478 GEN
 2479 lcmii(GEN x, GEN y)
 2480 {
 2481   pari_sp av;
 2482   GEN a, b;
 2483   if (!signe(x) || !signe(y)) return gen_0;
 2484   av = avma; a = gcdii(x,y);
 2485   if (absequalii(a,y)) { set_avma(av); return absi(x); }
 2486   if (!equali1(a)) y = diviiexact(y,a);
 2487   b = mulii(x,y); setabssign(b); return gerepileuptoint(av, b);
 2488 }
 2489 
 2490 /* given x in assume 0 < x < N; return u in (Z/NZ)^* such that u x = gcd(x,N) (mod N);
 2491  * set *pd = gcd(x,N) */
 2492 GEN
 2493 Fp_invgen(GEN x, GEN N, GEN *pd)
 2494 {
 2495   GEN d, d0, e, v;
 2496   if (lgefint(N) == 3)
 2497   {
 2498     ulong dd, NN = N[2], xx = umodiu(x,NN);
 2499     if (!xx) { *pd = N; return gen_0; }
 2500     xx = Fl_invgen(xx, NN, &dd);
 2501     *pd = utoi(dd); return utoi(xx);
 2502   }
 2503   *pd = d = bezout(x, N, &v, NULL);
 2504   if (equali1(d)) return v;
 2505   /* vx = gcd(x,N) (mod N), v coprime to N/d but need not be coprime to N */
 2506   e = diviiexact(N,d);
 2507   d0 = Z_ppo(d, e); /* d = d0 d1, d0 coprime to N/d, rad(d1) | N/d */
 2508   if (equali1(d0)) return v;
 2509   if (!equalii(d,d0)) e = lcmii(e, diviiexact(d,d0));
 2510   return Z_chinese_coprime(v, gen_1, e, d0, mulii(e,d0));
 2511 }
 2512 
 2513 /*********************************************************************/
 2514 /**                                                                 **/
 2515 /**                      CHINESE REMAINDERS                         **/
 2516 /**                                                                 **/
 2517 /*********************************************************************/
 2518 
 2519 /* Chinese Remainder Theorem.  x and y must have the same type (integermod,
 2520  * polymod, or polynomial/vector/matrix recursively constructed with these
 2521  * as coefficients). Creates (with the same type) a z in the same residue
 2522  * class as x and the same residue class as y, if it is possible.
 2523  *
 2524  * We also allow (during recursion) two identical objects even if they are
 2525  * not integermod or polymod. For example:
 2526  *
 2527  * ? x = [1, Mod(5, 11), Mod(X + Mod(2, 7), X^2 + 1)];
 2528  * ? y = [1, Mod(7, 17), Mod(X + Mod(0, 3), X^2 + 1)];
 2529  * ? chinese(x, y)
 2530  * %3 = [1, Mod(16, 187), Mod(X + mod(9, 21), X^2 + 1)] */
 2531 
 2532 static GEN
 2533 gen_chinese(GEN x, GEN(*f)(GEN,GEN))
 2534 {
 2535   GEN z = gassoc_proto(f,x,NULL);
 2536   if (z == gen_1) retmkintmod(gen_0,gen_1);
 2537   return z;
 2538 }
 2539 
 2540 /* x t_INTMOD, y t_POLMOD; promote x to t_POLMOD mod Pol(x.mod) then
 2541  * call chinese: makes Mod(0,1) a better "neutral" element */
 2542 static GEN
 2543 chinese_intpol(GEN x,GEN y)
 2544 {
 2545   pari_sp av = avma;
 2546   GEN z = mkpolmod(gel(x,2), scalarpol_shallow(gel(x,1), varn(gel(y,1))));
 2547   return gerepileupto(av, chinese(z, y));
 2548 }
 2549 
 2550 GEN
 2551 chinese1(GEN x) { return gen_chinese(x,chinese); }
 2552 
 2553 GEN
 2554 chinese(GEN x, GEN y)
 2555 {
 2556   pari_sp av;
 2557   long tx = typ(x), ty;
 2558   GEN z,p1,p2,d,u,v;
 2559 
 2560   if (!y) return chinese1(x);
 2561   if (gidentical(x,y)) return gcopy(x);
 2562   ty = typ(y);
 2563   if (tx == ty) switch(tx)
 2564   {
 2565     case t_POLMOD:
 2566     {
 2567       GEN A = gel(x,1), B = gel(y,1);
 2568       GEN a = gel(x,2), b = gel(y,2);
 2569       if (varn(A)!=varn(B)) pari_err_VAR("chinese",A,B);
 2570       if (RgX_equal(A,B)) retmkpolmod(chinese(a,b), gcopy(A)); /*same modulus*/
 2571       av = avma;
 2572       d = RgX_extgcd(A,B,&u,&v);
 2573       p2 = gsub(b, a);
 2574       if (!gequal0(gmod(p2, d))) break;
 2575       p1 = gdiv(A,d);
 2576       p2 = gadd(a, gmul(gmul(u,p1), p2));
 2577 
 2578       z = cgetg(3, t_POLMOD);
 2579       gel(z,1) = gmul(p1,B);
 2580       gel(z,2) = gmod(p2,gel(z,1));
 2581       return gerepileupto(av, z);
 2582     }
 2583     case t_INTMOD:
 2584     {
 2585       GEN A = gel(x,1), B = gel(y,1);
 2586       GEN a = gel(x,2), b = gel(y,2), c, d, C, U;
 2587       z = cgetg(3,t_INTMOD);
 2588       Z_chinese_pre(A, B, &C, &U, &d);
 2589       c = Z_chinese_post(a, b, C, U, d);
 2590       if (!c) pari_err_OP("chinese", x,y);
 2591       set_avma((pari_sp)z);
 2592       gel(z,1) = icopy(C);
 2593       gel(z,2) = icopy(c); return z;
 2594     }
 2595     case t_POL:
 2596     {
 2597       long i, lx = lg(x), ly = lg(y);
 2598       if (varn(x) != varn(y)) break;
 2599       if (lx < ly) { swap(x,y); lswap(lx,ly); }
 2600       z = cgetg(lx, t_POL); z[1] = x[1];
 2601       for (i=2; i<ly; i++) gel(z,i) = chinese(gel(x,i),gel(y,i));
 2602       for (   ; i<lx; i++) gel(z,i) = gcopy(gel(x,i));
 2603       return z;
 2604     }
 2605 
 2606     case t_VEC: case t_COL: case t_MAT:
 2607     {
 2608       long i, lx;
 2609       z = cgetg_copy(x, &lx); if (lx!=lg(y)) break;
 2610       for (i=1; i<lx; i++) gel(z,i) = chinese(gel(x,i),gel(y,i));
 2611       return z;
 2612     }
 2613   }
 2614   if (tx == t_POLMOD && ty == t_INTMOD) return chinese_intpol(y,x);
 2615   if (ty == t_POLMOD && tx == t_INTMOD) return chinese_intpol(x,y);
 2616   pari_err_OP("chinese",x,y);
 2617   return NULL; /* LCOV_EXCL_LINE */
 2618 }
 2619 
 2620 /* init chinese(Mod(.,A), Mod(.,B)) */
 2621 void
 2622 Z_chinese_pre(GEN A, GEN B, GEN *pC, GEN *pU, GEN *pd)
 2623 {
 2624   GEN u, d = bezout(A,B,&u,NULL); /* U = u(A/d), u(A/d) + v(B/d) = 1 */
 2625   GEN t = diviiexact(A,d);
 2626   *pU = mulii(u, t);
 2627   *pC = mulii(t, B);
 2628   if (pd) *pd = d;
 2629 }
 2630 /* Assume C = lcm(A, B), U = 0 mod (A/d), U = 1 mod (B/d), a = b mod d,
 2631  * where d = gcd(A,B) or NULL, return x = a (mod A), b (mod B).
 2632  * If d not NULL, check whether a = b mod d. */
 2633 GEN
 2634 Z_chinese_post(GEN a, GEN b, GEN C, GEN U, GEN d)
 2635 {
 2636   GEN b_a;
 2637   if (!signe(a))
 2638   {
 2639     if (d && !dvdii(b, d)) return NULL;
 2640     return Fp_mul(b, U, C);
 2641   }
 2642   b_a = subii(b,a);
 2643   if (d && !dvdii(b_a, d)) return NULL;
 2644   return modii(addii(a, mulii(U, b_a)), C);
 2645 }
 2646 static ulong
 2647 u_chinese_post(ulong a, ulong b, ulong C, ulong U)
 2648 {
 2649   if (!a) return Fl_mul(b, U, C);
 2650   return Fl_add(a, Fl_mul(U, Fl_sub(b,a,C), C), C);
 2651 }
 2652 
 2653 GEN
 2654 Z_chinese(GEN a, GEN b, GEN A, GEN B)
 2655 {
 2656   pari_sp av = avma;
 2657   GEN C, U; Z_chinese_pre(A, B, &C, &U, NULL);
 2658   return gerepileuptoint(av, Z_chinese_post(a,b, C, U, NULL));
 2659 }
 2660 GEN
 2661 Z_chinese_all(GEN a, GEN b, GEN A, GEN B, GEN *pC)
 2662 {
 2663   GEN U; Z_chinese_pre(A, B, pC, &U, NULL);
 2664   return Z_chinese_post(a,b, *pC, U, NULL);
 2665 }
 2666 
 2667 /* return lift(chinese(a mod A, b mod B))
 2668  * assume(A,B)=1, a,b,A,B integers and C = A*B */
 2669 GEN
 2670 Z_chinese_coprime(GEN a, GEN b, GEN A, GEN B, GEN C)
 2671 {
 2672   pari_sp av = avma;
 2673   GEN U = mulii(Fp_inv(A,B), A);
 2674   return gerepileuptoint(av, Z_chinese_post(a,b,C,U, NULL));
 2675 }
 2676 ulong
 2677 u_chinese_coprime(ulong a, ulong b, ulong A, ulong B, ulong C)
 2678 { return u_chinese_post(a,b,C, A * Fl_inv(A % B,B)); }
 2679 
 2680 /* chinese1 for coprime moduli in Z */
 2681 static GEN
 2682 chinese1_coprime_Z_aux(GEN x, GEN y)
 2683 {
 2684   GEN z = cgetg(3, t_INTMOD);
 2685   GEN A = gel(x,1), a = gel(x, 2);
 2686   GEN B = gel(y,1), b = gel(y, 2), C = mulii(A,B);
 2687   pari_sp av = avma;
 2688   GEN U = mulii(Fp_inv(A,B), A);
 2689   gel(z,2) = gerepileuptoint(av, Z_chinese_post(a,b,C,U, NULL));
 2690   gel(z,1) = C; return z;
 2691 }
 2692 GEN
 2693 chinese1_coprime_Z(GEN x) {return gen_chinese(x,chinese1_coprime_Z_aux);}
 2694 
 2695 /*********************************************************************/
 2696 /**                                                                 **/
 2697 /**                    MODULAR EXPONENTIATION                       **/
 2698 /**                                                                 **/
 2699 /*********************************************************************/
 2700 
 2701 /* xa, ya = t_VECSMALL */
 2702 GEN
 2703 ZV_producttree(GEN xa)
 2704 {
 2705   long n = lg(xa)-1;
 2706   long m = n==1 ? 1: expu(n-1)+1;
 2707   GEN T = cgetg(m+1, t_VEC), t;
 2708   long i, j, k;
 2709   t = cgetg(((n+1)>>1)+1, t_VEC);
 2710   if (typ(xa)==t_VECSMALL)
 2711   {
 2712     for (j=1, k=1; k<n; j++, k+=2)
 2713       gel(t, j) = muluu(xa[k], xa[k+1]);
 2714     if (k==n) gel(t, j) = utoi(xa[k]);
 2715   } else {
 2716     for (j=1, k=1; k<n; j++, k+=2)
 2717       gel(t, j) = mulii(gel(xa,k), gel(xa,k+1));
 2718     if (k==n) gel(t, j) = icopy(gel(xa,k));
 2719   }
 2720   gel(T,1) = t;
 2721   for (i=2; i<=m; i++)
 2722   {
 2723     GEN u = gel(T, i-1);
 2724     long n = lg(u)-1;
 2725     t = cgetg(((n+1)>>1)+1, t_VEC);
 2726     for (j=1, k=1; k<n; j++, k+=2)
 2727       gel(t, j) = mulii(gel(u, k), gel(u, k+1));
 2728     if (k==n) gel(t, j) = gel(u, k);
 2729     gel(T, i) = t;
 2730   }
 2731   return T;
 2732 }
 2733 
 2734 /* return [A mod P[i], i=1..#P], T = ZV_producttree(P) */
 2735 GEN
 2736 Z_ZV_mod_tree(GEN A, GEN P, GEN T)
 2737 {
 2738   long i,j,k;
 2739   long m = lg(T)-1, n = lg(P)-1;
 2740   GEN t;
 2741   GEN Tp = cgetg(m+1, t_VEC);
 2742   gel(Tp, m) = mkvec(A);
 2743   for (i=m-1; i>=1; i--)
 2744   {
 2745     GEN u = gel(T, i);
 2746     GEN v = gel(Tp, i+1);
 2747     long n = lg(u)-1;
 2748     t = cgetg(n+1, t_VEC);
 2749     for (j=1, k=1; k<n; j++, k+=2)
 2750     {
 2751       gel(t, k)   = modii(gel(v, j), gel(u, k));
 2752       gel(t, k+1) = modii(gel(v, j), gel(u, k+1));
 2753     }
 2754     if (k==n) gel(t, k) = gel(v, j);
 2755     gel(Tp, i) = t;
 2756   }
 2757   {
 2758     GEN u = gel(T, i+1);
 2759     GEN v = gel(Tp, i+1);
 2760     long l = lg(u)-1;
 2761     if (typ(P)==t_VECSMALL)
 2762     {
 2763       GEN R = cgetg(n+1, t_VECSMALL);
 2764       for (j=1, k=1; j<=l; j++, k+=2)
 2765       {
 2766         uel(R,k) = umodiu(gel(v, j), P[k]);
 2767         if (k < n)
 2768           uel(R,k+1) = umodiu(gel(v, j), P[k+1]);
 2769       }
 2770       return R;
 2771     }
 2772     else
 2773     {
 2774       GEN R = cgetg(n+1, t_VEC);
 2775       for (j=1, k=1; j<=l; j++, k+=2)
 2776       {
 2777         gel(R,k) = modii(gel(v, j), gel(P,k));
 2778         if (k < n)
 2779           gel(R,k+1) = modii(gel(v, j), gel(P,k+1));
 2780       }
 2781       return R;
 2782     }
 2783   }
 2784 }
 2785 
 2786 /* T = ZV_producttree(P), R = ZV_chinesetree(P,T) */
 2787 GEN
 2788 ZV_chinese_tree(GEN A, GEN P, GEN T, GEN R)
 2789 {
 2790   long m = lg(T)-1, n = lg(A)-1;
 2791   long i,j,k;
 2792   GEN Tp = cgetg(m+1, t_VEC);
 2793   GEN M = gel(T, 1);
 2794   GEN t = cgetg(lg(M), t_VEC);
 2795   if (typ(P)==t_VECSMALL)
 2796   {
 2797     for (j=1, k=1; k<n; j++, k+=2)
 2798     {
 2799       pari_sp av = avma;
 2800       GEN a = mului(A[k], gel(R,k)), b = mului(A[k+1], gel(R,k+1));
 2801       GEN tj = modii(addii(mului(P[k],b), mului(P[k+1],a)), gel(M,j));
 2802       gel(t, j) = gerepileuptoint(av, tj);
 2803     }
 2804     if (k==n) gel(t, j) = modii(mului(A[k], gel(R,k)), gel(M, j));
 2805   } else
 2806   {
 2807     for (j=1, k=1; k<n; j++, k+=2)
 2808     {
 2809       pari_sp av = avma;
 2810       GEN a = mulii(gel(A,k), gel(R,k)), b = mulii(gel(A,k+1), gel(R,k+1));
 2811       GEN tj = modii(addii(mulii(gel(P,k),b), mulii(gel(P,k+1),a)), gel(M,j));
 2812       gel(t, j) = gerepileuptoint(av, tj);
 2813     }
 2814     if (k==n) gel(t, j) = modii(mulii(gel(A,k), gel(R,k)), gel(M, j));
 2815   }
 2816   gel(Tp, 1) = t;
 2817   for (i=2; i<=m; i++)
 2818   {
 2819     GEN u = gel(T, i-1), M = gel(T, i);
 2820     GEN t = cgetg(lg(M), t_VEC);
 2821     GEN v = gel(Tp, i-1);
 2822     long n = lg(v)-1;
 2823     for (j=1, k=1; k<n; j++, k+=2)
 2824     {
 2825       pari_sp av = avma;
 2826       gel(t, j) = gerepileuptoint(av, modii(addii(mulii(gel(u, k), gel(v, k+1)),
 2827             mulii(gel(u, k+1), gel(v, k))), gel(M, j)));
 2828     }
 2829     if (k==n) gel(t, j) = gel(v, k);
 2830     gel(Tp, i) = t;
 2831   }
 2832   return gmael(Tp,m,1);
 2833 }
 2834 
 2835 static GEN
 2836 ncV_polint_center_tree(GEN vA, GEN P, GEN T, GEN R, GEN m2)
 2837 {
 2838   long i, l = lg(gel(vA,1)), n = lg(P);
 2839   GEN mod = gmael(T, lg(T)-1, 1), V = cgetg(l, t_COL);
 2840   for (i=1; i < l; i++)
 2841   {
 2842     pari_sp av = avma;
 2843     GEN c, A = cgetg(n, typ(P));
 2844     long j;
 2845     for (j=1; j < n; j++) A[j] = mael(vA,j,i);
 2846     c = Fp_center(ZV_chinese_tree(A, P, T, R), mod, m2);
 2847     gel(V,i) = gerepileuptoint(av, c);
 2848   }
 2849   return V;
 2850 }
 2851 
 2852 static GEN
 2853 nxV_polint_center_tree(GEN vA, GEN P, GEN T, GEN R, GEN m2)
 2854 {
 2855   long i, j, l, n = lg(P);
 2856   GEN mod = gmael(T, lg(T)-1, 1), V, w;
 2857   w = cgetg(n, t_VECSMALL);
 2858   for(j=1; j<n; j++) w[j] = lg(gel(vA,j));
 2859   l = vecsmall_max(w);
 2860   V = cgetg(l, t_POL);
 2861   V[1] = mael(vA,1,1);
 2862   for (i=2; i < l; i++)
 2863   {
 2864     pari_sp av = avma;
 2865     GEN c, A = cgetg(n, typ(P));
 2866     if (typ(P)==t_VECSMALL)
 2867       for (j=1; j < n; j++) A[j] = i < w[j] ? mael(vA,j,i): 0;
 2868     else
 2869       for (j=1; j < n; j++) gel(A,j) = i < w[j] ? gmael(vA,j,i): gen_0;
 2870     c = Fp_center(ZV_chinese_tree(A, P, T, R), mod, m2);
 2871     gel(V,i) = gerepileuptoint(av, c);
 2872   }
 2873   return ZX_renormalize(V, l);
 2874 }
 2875 
 2876 static GEN
 2877 nxCV_polint_center_tree(GEN vA, GEN P, GEN T, GEN R, GEN m2)
 2878 {
 2879   long i, j, l = lg(gel(vA,1)), n = lg(P);
 2880   GEN A = cgetg(n, t_VEC);
 2881   GEN V = cgetg(l, t_COL);
 2882   for (i=1; i < l; i++)
 2883   {
 2884     for (j=1; j < n; j++) gel(A,j) = gmael(vA,j,i);
 2885     gel(V,i) = nxV_polint_center_tree(A, P, T, R, m2);
 2886   }
 2887   return V;
 2888 }
 2889 
 2890 static GEN
 2891 polint_chinese(GEN worker, GEN mA, GEN P)
 2892 {
 2893   long cnt, pending, n, i, j, l = lg(gel(mA,1));
 2894   struct pari_mt pt;
 2895   GEN done, va, M, A;
 2896   pari_timer ti;
 2897 
 2898   if (l == 1) return cgetg(1, t_MAT);
 2899   cnt = pending = 0;
 2900   n = lg(P);
 2901   A = cgetg(n, t_VEC);
 2902   va = mkvec(A);
 2903   M = cgetg(l, t_MAT);
 2904   if (DEBUGLEVEL>4) timer_start(&ti);
 2905   if (DEBUGLEVEL>5) err_printf("Start parallel Chinese remainder: ");
 2906   mt_queue_start_lim(&pt, worker, l-1);
 2907   for (i=1; i<l || pending; i++)
 2908   {
 2909     long workid;
 2910     for(j=1; j < n; j++) gel(A,j) = gmael(mA,j,i);
 2911     mt_queue_submit(&pt, i, i<l? va: NULL);
 2912     done = mt_queue_get(&pt, &workid, &pending);
 2913     if (done)
 2914     {
 2915       gel(M,workid) = done;
 2916       if (DEBUGLEVEL>5) err_printf("%ld%% ",(++cnt)*100/(l-1));
 2917     }
 2918   }
 2919   if (DEBUGLEVEL>5) err_printf("\n");
 2920   if (DEBUGLEVEL>4) timer_printf(&ti, "nmV_chinese_center");
 2921   mt_queue_end(&pt);
 2922   return M;
 2923 }
 2924 
 2925 GEN
 2926 nxMV_polint_center_tree_worker(GEN vA, GEN T, GEN R, GEN P, GEN m2)
 2927 {
 2928   return nxCV_polint_center_tree(vA, P, T, R, m2);
 2929 }
 2930 
 2931 static GEN
 2932 nxMV_polint_center_tree_seq(GEN vA, GEN P, GEN T, GEN R, GEN m2)
 2933 {
 2934   long i, j, l = lg(gel(vA,1)), n = lg(P);
 2935   GEN A = cgetg(n, t_VEC);
 2936   GEN V = cgetg(l, t_MAT);
 2937   for (i=1; i < l; i++)
 2938   {
 2939     for (j=1; j < n; j++) gel(A,j) = gmael(vA,j,i);
 2940     gel(V,i) = nxCV_polint_center_tree(A, P, T, R, m2);
 2941   }
 2942   return V;
 2943 }
 2944 
 2945 static GEN
 2946 nxMV_polint_center_tree(GEN mA, GEN P, GEN T, GEN R, GEN m2)
 2947 {
 2948   GEN worker = snm_closure(is_entry("_nxMV_polint_worker"), mkvec4(T, R, P, m2));
 2949   return polint_chinese(worker, mA, P);
 2950 }
 2951 
 2952 static GEN
 2953 nmV_polint_center_tree_seq(GEN vA, GEN P, GEN T, GEN R, GEN m2)
 2954 {
 2955   long i, j, l = lg(gel(vA,1)), n = lg(P);
 2956   GEN A = cgetg(n, t_VEC);
 2957   GEN V = cgetg(l, t_MAT);
 2958   for (i=1; i < l; i++)
 2959   {
 2960     for (j=1; j < n; j++) gel(A,j) = gmael(vA,j,i);
 2961     gel(V,i) = ncV_polint_center_tree(A, P, T, R, m2);
 2962   }
 2963   return V;
 2964 }
 2965 
 2966 GEN
 2967 nmV_polint_center_tree_worker(GEN vA, GEN T, GEN R, GEN P, GEN m2)
 2968 {
 2969   return ncV_polint_center_tree(vA, P, T, R, m2);
 2970 }
 2971 
 2972 static GEN
 2973 nmV_polint_center_tree(GEN mA, GEN P, GEN T, GEN R, GEN m2)
 2974 {
 2975   GEN worker = snm_closure(is_entry("_polint_worker"), mkvec4(T, R, P, m2));
 2976   return polint_chinese(worker, mA, P);
 2977 }
 2978 
 2979 /* return [A mod P[i], i=1..#P] */
 2980 GEN
 2981 Z_ZV_mod(GEN A, GEN P)
 2982 {
 2983   pari_sp av = avma;
 2984   return gerepilecopy(av, Z_ZV_mod_tree(A, P, ZV_producttree(P)));
 2985 }
 2986 /* P a t_VECSMALL */
 2987 GEN
 2988 Z_nv_mod(GEN A, GEN P)
 2989 {
 2990   pari_sp av = avma;
 2991   return gerepileuptoleaf(av, Z_ZV_mod_tree(A, P, ZV_producttree(P)));
 2992 }
 2993 /* B a ZX, T = ZV_producttree(P) */
 2994 GEN
 2995 ZX_nv_mod_tree(GEN B, GEN A, GEN T)
 2996 {
 2997   pari_sp av;
 2998   long i, j, l = lg(B), n = lg(A)-1;
 2999   GEN V = cgetg(n+1, t_VEC);
 3000   for (j=1; j <= n; j++)
 3001   {
 3002     gel(V, j) = cgetg(l, t_VECSMALL);
 3003     mael(V, j, 1) = B[1]&VARNBITS;
 3004   }
 3005   av = avma;
 3006   for (i=2; i < l; i++)
 3007   {
 3008     GEN v = Z_ZV_mod_tree(gel(B, i), A, T);
 3009     for (j=1; j <= n; j++)
 3010       mael(V, j, i) = v[j];
 3011     set_avma(av);
 3012   }
 3013   for (j=1; j <= n; j++)
 3014     (void) Flx_renormalize(gel(V, j), l);
 3015   return V;
 3016 }
 3017 
 3018 static GEN
 3019 to_ZX(GEN a, long v) { return typ(a)==t_INT? scalarpol(a,v): a; }
 3020 
 3021 GEN
 3022 ZXX_nv_mod_tree(GEN P, GEN xa, GEN T, long w)
 3023 {
 3024   pari_sp av = avma;
 3025   long i, j, l = lg(P), n = lg(xa)-1, vP = varn(P);
 3026   GEN V = cgetg(n+1, t_VEC);
 3027   for (j=1; j <= n; j++)
 3028   {
 3029     gel(V, j) = cgetg(l, t_POL);
 3030     mael(V, j, 1) = vP;
 3031   }
 3032   for (i=2; i < l; i++)
 3033   {
 3034     GEN v = ZX_nv_mod_tree(to_ZX(gel(P, i), w), xa, T);
 3035     for (j=1; j <= n; j++)
 3036       gmael(V, j, i) = gel(v,j);
 3037   }
 3038   for (j=1; j <= n; j++)
 3039     (void) FlxX_renormalize(gel(V, j), l);
 3040   return gerepilecopy(av, V);
 3041 }
 3042 
 3043 GEN
 3044 ZXC_nv_mod_tree(GEN C, GEN xa, GEN T, long w)
 3045 {
 3046   pari_sp av = avma;
 3047   long i, j, l = lg(C), n = lg(xa)-1;
 3048   GEN V = cgetg(n+1, t_VEC);
 3049   for (j = 1; j <= n; j++)
 3050     gel(V, j) = cgetg(l, t_COL);
 3051   for (i = 1; i < l; i++)
 3052   {
 3053     GEN v = ZX_nv_mod_tree(to_ZX(gel(C, i), w), xa, T);
 3054     for (j = 1; j <= n; j++)
 3055       gmael(V, j, i) = gel(v,j);
 3056   }
 3057   return gerepilecopy(av, V);
 3058 }
 3059 
 3060 GEN
 3061 ZXM_nv_mod_tree(GEN M, GEN xa, GEN T, long w)
 3062 {
 3063   pari_sp av = avma;
 3064   long i, j, l = lg(M), n = lg(xa)-1;
 3065   GEN V = cgetg(n+1, t_VEC);
 3066   for (j=1; j <= n; j++)
 3067     gel(V, j) = cgetg(l, t_MAT);
 3068   for (i=1; i < l; i++)
 3069   {
 3070     GEN v = ZXC_nv_mod_tree(gel(M, i), xa, T, w);
 3071     for (j=1; j <= n; j++)
 3072       gmael(V, j, i) = gel(v,j);
 3073   }
 3074   return gerepilecopy(av, V);
 3075 }
 3076 
 3077 GEN
 3078 ZV_nv_mod_tree(GEN B, GEN A, GEN T)
 3079 {
 3080   pari_sp av;
 3081   long i, j, l = lg(B), n = lg(A)-1;
 3082   GEN V = cgetg(n+1, t_VEC);
 3083   for (j=1; j <= n; j++)
 3084     gel(V, j) = cgetg(l, t_VECSMALL);
 3085   av = avma;
 3086   for (i=1; i < l; i++)
 3087   {
 3088     GEN v = Z_ZV_mod_tree(gel(B, i), A, T);
 3089     for (j=1; j <= n; j++)
 3090       mael(V, j, i) = v[j];
 3091     set_avma(av);
 3092   }
 3093   return V;
 3094 }
 3095 
 3096 GEN
 3097 ZM_nv_mod_tree(GEN M, GEN xa, GEN T)
 3098 {
 3099   pari_sp av = avma;
 3100   long i, j, l = lg(M), n = lg(xa)-1;
 3101   GEN V = cgetg(n+1, t_VEC);
 3102   for (j=1; j <= n; j++)
 3103     gel(V, j) = cgetg(l, t_MAT);
 3104   for (i=1; i < l; i++)
 3105   {
 3106     GEN v = ZV_nv_mod_tree(gel(M, i), xa, T);
 3107     for (j=1; j <= n; j++)
 3108       gmael(V, j, i) = gel(v,j);
 3109   }
 3110   return gerepilecopy(av, V);
 3111 }
 3112 
 3113 static GEN
 3114 ZV_sqr(GEN z)
 3115 {
 3116   long i,l = lg(z);
 3117   GEN x = cgetg(l, t_VEC);
 3118   if (typ(z)==t_VECSMALL)
 3119     for (i=1; i<l; i++) gel(x,i) = sqru(z[i]);
 3120   else
 3121     for (i=1; i<l; i++) gel(x,i) = sqri(gel(z,i));
 3122   return x;
 3123 }
 3124 
 3125 static GEN
 3126 ZT_sqr(GEN x)
 3127 {
 3128   if (typ(x) == t_INT)
 3129     return sqri(x);
 3130   pari_APPLY_type(t_VEC, ZT_sqr(gel(x,i)))
 3131 }
 3132 
 3133 static GEN
 3134 ZV_invdivexact(GEN y, GEN x)
 3135 {
 3136   long i, l = lg(y);
 3137   GEN z = cgetg(l,t_VEC);
 3138   if (typ(x)==t_VECSMALL)
 3139     for (i=1; i<l; i++)
 3140     {
 3141       pari_sp av = avma;
 3142       ulong a = Fl_inv(umodiu(diviuexact(gel(y,i),x[i]), x[i]), x[i]);
 3143       set_avma(av);
 3144       gel(z,i) = utoipos(a);
 3145     }
 3146   else
 3147     for (i=1; i<l; i++)
 3148       gel(z,i) = Fp_inv(diviiexact(gel(y,i), gel(x,i)), gel(x,i));
 3149   return z;
 3150 }
 3151 
 3152 /* P t_VECSMALL or t_VEC of t_INT  */
 3153 GEN
 3154 ZV_chinesetree(GEN P, GEN T)
 3155 {
 3156   GEN T2 = ZT_sqr(T), P2 = ZV_sqr(P);
 3157   GEN mod = gmael(T,lg(T)-1,1);
 3158   return ZV_invdivexact(Z_ZV_mod_tree(mod, P2, T2), P);
 3159 }
 3160 
 3161 static GEN
 3162 gc_chinese(pari_sp av, GEN T, GEN a, GEN *pt_mod)
 3163 {
 3164   if (!pt_mod)
 3165     return gerepileupto(av, a);
 3166   else
 3167   {
 3168     GEN mod = gmael(T, lg(T)-1, 1);
 3169     gerepileall(av, 2, &a, &mod);
 3170     *pt_mod = mod;
 3171     return a;
 3172   }
 3173 }
 3174 
 3175 GEN
 3176 ZV_chinese_center(GEN A, GEN P, GEN *pt_mod)
 3177 {
 3178   pari_sp av = avma;
 3179   GEN T = ZV_producttree(P);
 3180   GEN R = ZV_chinesetree(P, T);
 3181   GEN a = ZV_chinese_tree(A, P, T, R);
 3182   GEN mod = gmael(T, lg(T)-1, 1);
 3183   GEN ca = Fp_center(a, mod, shifti(mod,-1));
 3184   return gc_chinese(av, T, ca, pt_mod);
 3185 }
 3186 
 3187 GEN
 3188 ZV_chinese(GEN A, GEN P, GEN *pt_mod)
 3189 {
 3190   pari_sp av = avma;
 3191   GEN T = ZV_producttree(P);
 3192   GEN R = ZV_chinesetree(P, T);
 3193   GEN a = ZV_chinese_tree(A, P, T, R);
 3194   return gc_chinese(av, T, a, pt_mod);
 3195 }
 3196 
 3197 GEN
 3198 nxV_chinese_center_tree(GEN A, GEN P, GEN T, GEN R)
 3199 {
 3200   pari_sp av = avma;
 3201   GEN m2 = shifti(gmael(T, lg(T)-1, 1), -1);
 3202   GEN a = nxV_polint_center_tree(A, P, T, R, m2);
 3203   return gerepileupto(av, a);
 3204 }
 3205 
 3206 GEN
 3207 nxV_chinese_center(GEN A, GEN P, GEN *pt_mod)
 3208 {
 3209   pari_sp av = avma;
 3210   GEN T = ZV_producttree(P);
 3211   GEN R = ZV_chinesetree(P, T);
 3212   GEN m2 = shifti(gmael(T, lg(T)-1, 1), -1);
 3213   GEN a = nxV_polint_center_tree(A, P, T, R, m2);
 3214   return gc_chinese(av, T, a, pt_mod);
 3215 }
 3216 
 3217 GEN
 3218 ncV_chinese_center(GEN A, GEN P, GEN *pt_mod)
 3219 {
 3220   pari_sp av = avma;
 3221   GEN T = ZV_producttree(P);
 3222   GEN R = ZV_chinesetree(P, T);
 3223   GEN m2 = shifti(gmael(T, lg(T)-1, 1), -1);
 3224   GEN a = ncV_polint_center_tree(A, P, T, R, m2);
 3225   return gc_chinese(av, T, a, pt_mod);
 3226 }
 3227 
 3228 GEN
 3229 ncV_chinese_center_tree(GEN A, GEN P, GEN T, GEN R)
 3230 {
 3231   pari_sp av = avma;
 3232   GEN m2 = shifti(gmael(T, lg(T)-1, 1), -1);
 3233   GEN a = ncV_polint_center_tree(A, P, T, R, m2);
 3234   return gerepileupto(av, a);
 3235 }
 3236 
 3237 GEN
 3238 nmV_chinese_center_tree(GEN A, GEN P, GEN T, GEN R)
 3239 {
 3240   pari_sp av = avma;
 3241   GEN m2 = shifti(gmael(T, lg(T)-1, 1), -1);
 3242   GEN a = nmV_polint_center_tree(A, P, T, R, m2);
 3243   return gerepileupto(av, a);
 3244 }
 3245 
 3246 GEN
 3247 nmV_chinese_center_tree_seq(GEN A, GEN P, GEN T, GEN R)
 3248 {
 3249   pari_sp av = avma;
 3250   GEN m2 = shifti(gmael(T, lg(T)-1, 1), -1);
 3251   GEN a = nmV_polint_center_tree_seq(A, P, T, R, m2);
 3252   return gerepileupto(av, a);
 3253 }
 3254 
 3255 GEN
 3256 nmV_chinese_center(GEN A, GEN P, GEN *pt_mod)
 3257 {
 3258   pari_sp av = avma;
 3259   GEN T = ZV_producttree(P);
 3260   GEN R = ZV_chinesetree(P, T);
 3261   GEN m2 = shifti(gmael(T, lg(T)-1, 1), -1);
 3262   GEN a = nmV_polint_center_tree(A, P, T, R, m2);
 3263   return gc_chinese(av, T, a, pt_mod);
 3264 }
 3265 
 3266 GEN
 3267 nxCV_chinese_center_tree(GEN A, GEN P, GEN T, GEN R)
 3268 {
 3269   pari_sp av = avma;
 3270   GEN m2 = shifti(gmael(T, lg(T)-1, 1), -1);
 3271   GEN a = nxCV_polint_center_tree(A, P, T, R, m2);
 3272   return gerepileupto(av, a);
 3273 }
 3274 
 3275 GEN
 3276 nxCV_chinese_center(GEN A, GEN P, GEN *pt_mod)
 3277 {
 3278   pari_sp av = avma;
 3279   GEN T = ZV_producttree(P);
 3280   GEN R = ZV_chinesetree(P, T);
 3281   GEN m2 = shifti(gmael(T, lg(T)-1, 1), -1);
 3282   GEN a = nxCV_polint_center_tree(A, P, T, R, m2);
 3283   return gc_chinese(av, T, a, pt_mod);
 3284 }
 3285 
 3286 GEN
 3287 nxMV_chinese_center_tree_seq(GEN A, GEN P, GEN T, GEN R)
 3288 {
 3289   pari_sp av = avma;
 3290   GEN m2 = shifti(gmael(T, lg(T)-1, 1), -1);
 3291   GEN a = nxMV_polint_center_tree_seq(A, P, T, R, m2);
 3292   return gerepileupto(av, a);
 3293 }
 3294 
 3295 GEN
 3296 nxMV_chinese_center(GEN A, GEN P, GEN *pt_mod)
 3297 {
 3298   pari_sp av = avma;
 3299   GEN T = ZV_producttree(P);
 3300   GEN R = ZV_chinesetree(P, T);
 3301   GEN m2 = shifti(gmael(T, lg(T)-1, 1), -1);
 3302   GEN a = nxMV_polint_center_tree(A, P, T, R, m2);
 3303   return gc_chinese(av, T, a, pt_mod);
 3304 }
 3305 
 3306 /**********************************************************************
 3307  **                                                                  **
 3308  **                    Powering  over (Z/NZ)^*, small N              **
 3309  **                                                                  **
 3310  **********************************************************************/
 3311 
 3312 /* 2^n mod p; assume n > 1 */
 3313 static ulong
 3314 Fl_2powu_pre(ulong n, ulong p, ulong pi)
 3315 {
 3316   ulong y = 2;
 3317   int j = 1+bfffo(n);
 3318   /* normalize, i.e set highest bit to 1 (we know n != 0) */
 3319   n<<=j; j = BITS_IN_LONG-j; /* first bit is now implicit */
 3320   for (; j; n<<=1,j--)
 3321   {
 3322     y = Fl_sqr_pre(y,p,pi);
 3323     if (n & HIGHBIT) y = Fl_double(y, p);
 3324   }
 3325   return y;
 3326 }
 3327 
 3328 /* 2^n mod p; assume n > 1 and !(p & HIGHMASK) */
 3329 static ulong
 3330 Fl_2powu(ulong n, ulong p)
 3331 {
 3332   ulong y = 2;
 3333   int j = 1+bfffo(n);
 3334   /* normalize, i.e set highest bit to 1 (we know n != 0) */
 3335   n<<=j; j = BITS_IN_LONG-j; /* first bit is now implicit */
 3336   for (; j; n<<=1,j--)
 3337   {
 3338     y = (y*y) % p;
 3339     if (n & HIGHBIT) y = Fl_double(y, p);
 3340   }
 3341   return y;
 3342 }
 3343 
 3344 ulong
 3345 Fl_powu_pre(ulong x, ulong n0, ulong p, ulong pi)
 3346 {
 3347   ulong y, z, n;
 3348   if (n0 <= 1)
 3349   { /* frequent special cases */
 3350     if (n0 == 1) return x;
 3351     if (n0 == 0) return 1;
 3352   }
 3353   if (x <= 2)
 3354   {
 3355     if (x == 2) return Fl_2powu_pre(n0, p, pi);
 3356     return x; /* 0 or 1 */
 3357   }
 3358   y = 1; z = x; n = n0;
 3359   for(;;)
 3360   {
 3361     if (n&1) y = Fl_mul_pre(y,z,p,pi);
 3362     n>>=1; if (!n) return y;
 3363     z = Fl_sqr_pre(z,p,pi);
 3364   }
 3365 }
 3366 
 3367 ulong
 3368 Fl_powu(ulong x, ulong n0, ulong p)
 3369 {
 3370   ulong y, z, n;
 3371   if (n0 <= 2)
 3372   { /* frequent special cases */
 3373     if (n0 == 2) return Fl_sqr(x,p);
 3374     if (n0 == 1) return x;
 3375     if (n0 == 0) return 1;
 3376   }
 3377   if (x <= 1) return x; /* 0 or 1 */
 3378   if (p & HIGHMASK)
 3379     return Fl_powu_pre(x, n0, p, get_Fl_red(p));
 3380   if (x == 2) return Fl_2powu(n0, p);
 3381   y = 1; z = x; n = n0;
 3382   for(;;)
 3383   {
 3384     if (n&1) y = (y*z) % p;
 3385     n>>=1; if (!n) return y;
 3386     z = (z*z) % p;
 3387   }
 3388 }
 3389 
 3390 /* Reduce data dependency to maximize internal parallelism */
 3391 GEN
 3392 Fl_powers_pre(ulong x, long n, ulong p, ulong pi)
 3393 {
 3394   long i, k;
 3395   GEN powers = cgetg(n + 2, t_VECSMALL);
 3396   powers[1] = 1; if (n == 0) return powers;
 3397   powers[2] = x;
 3398   for (i = 3, k=2; i <= n; i+=2, k++)
 3399   {
 3400     powers[i] = Fl_sqr_pre(powers[k], p, pi);
 3401     powers[i+1] = Fl_mul_pre(powers[k], powers[k+1], p, pi);
 3402   }
 3403   if (i==n+1)
 3404     powers[i] = Fl_sqr_pre(powers[k], p, pi);
 3405   return powers;
 3406 }
 3407 
 3408 GEN
 3409 Fl_powers(ulong x, long n, ulong p)
 3410 {
 3411   return Fl_powers_pre(x, n, p, get_Fl_red(p));
 3412 }
 3413 
 3414 /**********************************************************************
 3415  **                                                                  **
 3416  **                    Powering  over (Z/NZ)^*, large N              **
 3417  **                                                                  **
 3418  **********************************************************************/
 3419 
 3420 static GEN
 3421 Fp_dblsqr(GEN x, GEN N)
 3422 {
 3423   GEN z = shifti(Fp_sqr(x, N), 1);
 3424   return cmpii(z, N) >= 0? subii(z, N): z;
 3425 }
 3426 
 3427 typedef struct muldata {
 3428   GEN (*sqr)(void * E, GEN x);
 3429   GEN (*mul)(void * E, GEN x, GEN y);
 3430   GEN (*mul2)(void * E, GEN x);
 3431 } muldata;
 3432 
 3433 /* modified Barrett reduction with one fold */
 3434 /* See Fast Modular Reduction, W. Hasenplaugh, G. Gaubatz, V. Gopal, ARITH 18 */
 3435 
 3436 static GEN
 3437 Fp_invmBarrett(GEN p, long s)
 3438 {
 3439   GEN R, Q = dvmdii(int2n(3*s),p,&R);
 3440   return mkvec2(Q,R);
 3441 }
 3442 
 3443 /* a <= (N-1)^2, 2^(2s-2) <= N < 2^(2s). Return 0 <= r < N such that
 3444  * a = r (mod N) */
 3445 static GEN
 3446 Fp_rem_mBarrett(GEN a, GEN B, long s, GEN N)
 3447 {
 3448   pari_sp av = avma;
 3449   GEN P = gel(B, 1), Q = gel(B, 2); /* 2^(3s) = P N + Q, 0 <= Q < N */
 3450   long t = expi(P)+1; /* 2^(t-1) <= P < 2^t */
 3451   GEN u = shifti(a, -3*s), v = remi2n(a, 3*s); /* a = 2^(3s)u + v */
 3452   GEN A = addii(v, mulii(Q,u)); /* 0 <= A < 2^(3s+1) */
 3453   GEN q = shifti(mulii(shifti(A, t-3*s), P), -t); /* A/N - 4 < q <= A/N */
 3454   GEN r = subii(A, mulii(q, N));
 3455   GEN sr= subii(r,N);     /* 0 <= r < 4*N */
 3456   if (signe(sr)<0) return gerepileuptoint(av, r);
 3457   r=sr; sr = subii(r,N);  /* 0 <= r < 3*N */
 3458   if (signe(sr)<0) return gerepileuptoint(av, r);
 3459   r=sr; sr = subii(r,N);  /* 0 <= r < 2*N */
 3460   return gerepileuptoint(av, signe(sr)>=0 ? sr:r);
 3461 }
 3462 
 3463 /* Montgomery reduction */
 3464 
 3465 INLINE ulong
 3466 init_montdata(GEN N) { return (ulong) -invmod2BIL(mod2BIL(N)); }
 3467 
 3468 struct montred
 3469 {
 3470   GEN N;
 3471   ulong inv;
 3472 };
 3473 
 3474 /* Montgomery reduction */
 3475 static GEN
 3476 _sqr_montred(void * E, GEN x)
 3477 {
 3478   struct montred * D = (struct montred *) E;
 3479   return red_montgomery(sqri(x), D->N, D->inv);
 3480 }
 3481 
 3482 /* Montgomery reduction */
 3483 static GEN
 3484 _mul_montred(void * E, GEN x, GEN y)
 3485 {
 3486   struct montred * D = (struct montred *) E;
 3487   return red_montgomery(mulii(x, y), D->N, D->inv);
 3488 }
 3489 
 3490 static GEN
 3491 _mul2_montred(void * E, GEN x)
 3492 {
 3493   struct montred * D = (struct montred *) E;
 3494   GEN z = shifti(_sqr_montred(E, x), 1);
 3495   long l = lgefint(D->N);
 3496   while (lgefint(z) > l) z = subii(z, D->N);
 3497   return z;
 3498 }
 3499 
 3500 static GEN
 3501 _sqr_remii(void* N, GEN x)
 3502 { return remii(sqri(x), (GEN) N); }
 3503 
 3504 static GEN
 3505 _mul_remii(void* N, GEN x, GEN y)
 3506 { return remii(mulii(x, y), (GEN) N); }
 3507 
 3508 static GEN
 3509 _mul2_remii(void* N, GEN x)
 3510 { return Fp_dblsqr(x, (GEN) N); }
 3511 
 3512 struct redbarrett
 3513 {
 3514   GEN iM, N;
 3515   long s;
 3516 };
 3517 
 3518 static GEN
 3519 _sqr_remiibar(void *E, GEN x)
 3520 {
 3521   struct redbarrett * D = (struct redbarrett *) E;
 3522   return Fp_rem_mBarrett(sqri(x), D->iM, D->s, D->N);
 3523 }
 3524 
 3525 static GEN
 3526 _mul_remiibar(void *E, GEN x, GEN y)
 3527 {
 3528   struct redbarrett * D = (struct redbarrett *) E;
 3529   return Fp_rem_mBarrett(mulii(x, y), D->iM, D->s, D->N);
 3530 }
 3531 
 3532 static GEN
 3533 _mul2_remiibar(void *E, GEN x)
 3534 {
 3535   struct redbarrett * D = (struct redbarrett *) E;
 3536   return Fp_dblsqr(x, D->N);
 3537 }
 3538 
 3539 static long
 3540 Fp_select_red(GEN *y, ulong k, GEN N, long lN, muldata *D, void **pt_E)
 3541 {
 3542   if (lN >= Fp_POW_BARRETT_LIMIT && (k==0 || ((double)k)*expi(*y) > 2 + expi(N)))
 3543   {
 3544     struct redbarrett * E = (struct redbarrett *) stack_malloc(sizeof(struct redbarrett));
 3545     D->sqr = &_sqr_remiibar;
 3546     D->mul = &_mul_remiibar;
 3547     D->mul2 = &_mul2_remiibar;
 3548     E->N = N;
 3549     E->s = 1+(expi(N)>>1);
 3550     E->iM = Fp_invmBarrett(N, E->s);
 3551     *pt_E = (void*) E;
 3552     return 0;
 3553   }
 3554   else if (mod2(N) && lN < Fp_POW_REDC_LIMIT)
 3555   {
 3556     struct montred * E = (struct montred *) stack_malloc(sizeof(struct montred));
 3557     *y = remii(shifti(*y, bit_accuracy(lN)), N);
 3558     D->sqr = &_sqr_montred;
 3559     D->mul = &_mul_montred;
 3560     D->mul2 = &_mul2_montred;
 3561     E->N = N;
 3562     E->inv = init_montdata(N);
 3563     *pt_E = (void*) E;
 3564     return 1;
 3565   }
 3566   else
 3567   {
 3568     D->sqr = &_sqr_remii;
 3569     D->mul = &_mul_remii;
 3570     D->mul2 = &_mul2_remii;
 3571     *pt_E = (void*) N;
 3572     return 0;
 3573   }
 3574 }
 3575 
 3576 GEN
 3577 Fp_powu(GEN A, ulong k, GEN N)
 3578 {
 3579   long lN = lgefint(N);
 3580   int base_is_2, use_montgomery;
 3581   muldata D;
 3582   void *E;
 3583   pari_sp av;
 3584 
 3585   if (lN == 3) {
 3586     ulong n = uel(N,2);
 3587     return utoi( Fl_powu(umodiu(A, n), k, n) );
 3588   }
 3589   if (k <= 2)
 3590   { /* frequent special cases */
 3591     if (k == 2) return Fp_sqr(A,N);
 3592     if (k == 1) return A;
 3593     if (k == 0) return gen_1;
 3594   }
 3595   av = avma; A = modii(A,N);
 3596   base_is_2 = 0;
 3597   if (lgefint(A) == 3) switch(A[2])
 3598   {
 3599     case 1: set_avma(av); return gen_1;
 3600     case 2:  base_is_2 = 1; break;
 3601   }
 3602 
 3603   /* TODO: Move this out of here and use for general modular computations */
 3604   use_montgomery = Fp_select_red(&A, k, N, lN, &D, &E);
 3605   if (base_is_2)
 3606     A = gen_powu_fold_i(A, k, E, D.sqr, D.mul2);
 3607   else
 3608     A = gen_powu_i(A, k, E, D.sqr, D.mul);
 3609   if (use_montgomery)
 3610   {
 3611     A = red_montgomery(A, N, ((struct montred *) E)->inv);
 3612     if (cmpii(A, N) >= 0) A = subii(A,N);
 3613   }
 3614   return gerepileuptoint(av, A);
 3615 }
 3616 
 3617 GEN
 3618 Fp_pows(GEN A, long k, GEN N)
 3619 {
 3620   if (lgefint(N) == 3) {
 3621     ulong n = N[2];
 3622     ulong a = umodiu(A, n);
 3623     if (k < 0) {
 3624       a = Fl_inv(a, n);
 3625       k = -k;
 3626     }
 3627     return utoi( Fl_powu(a, (ulong)k, n) );
 3628   }
 3629   if (k < 0) { A = Fp_inv(A, N); k = -k; };
 3630   return Fp_powu(A, (ulong)k, N);
 3631 }
 3632 
 3633 /* A^K mod N */
 3634 GEN
 3635 Fp_pow(GEN A, GEN K, GEN N)
 3636 {
 3637   pari_sp av;
 3638   long s, lN = lgefint(N), sA, sy;
 3639   int base_is_2, use_montgomery;
 3640   GEN y;
 3641   muldata D;
 3642   void *E;
 3643 
 3644   s = signe(K);
 3645   if (!s) return dvdii(A,N)? gen_0: gen_1;
 3646   if (lN == 3 && lgefint(K) == 3)
 3647   {
 3648     ulong n = N[2], a = umodiu(A, n);
 3649     if (s < 0) a = Fl_inv(a, n);
 3650     if (a <= 1) return utoi(a); /* 0 or 1 */
 3651     return utoi(Fl_powu(a, uel(K,2), n));
 3652   }
 3653 
 3654   av = avma;
 3655   if (s < 0) y = Fp_inv(A,N);
 3656   else
 3657   {
 3658     y = modii(A,N);
 3659     if (!signe(y)) { set_avma(av); return gen_0; }
 3660   }
 3661   if (lgefint(K) == 3) return gerepileuptoint(av, Fp_powu(y, K[2], N));
 3662 
 3663   base_is_2 = 0;
 3664   sy = abscmpii(y, shifti(N,-1)) > 0;
 3665   if (sy) y = subii(N,y);
 3666   sA = sy && mod2(K);
 3667   if (lgefint(y) == 3) switch(y[2])
 3668   {
 3669     case 1:  set_avma(av); return sA ? subis(N,1): gen_1;
 3670     case 2:  base_is_2 = 1; break;
 3671   }
 3672 
 3673   /* TODO: Move this out of here and use for general modular computations */
 3674   use_montgomery = Fp_select_red(&y, 0UL, N, lN, &D, &E);
 3675   if (base_is_2)
 3676     y = gen_pow_fold_i(y, K, E, D.sqr, D.mul2);
 3677   else
 3678     y = gen_pow_i(y, K, E, D.sqr, D.mul);
 3679   if (use_montgomery)
 3680   {
 3681     y = red_montgomery(y, N, ((struct montred *) E)->inv);
 3682     if (cmpii(y,N) >= 0) y = subii(y,N);
 3683   }
 3684   if (sA) y = subii(N, y);
 3685   return gerepileuptoint(av,y);
 3686 }
 3687 
 3688 static GEN
 3689 _Fp_mul(void *E, GEN x, GEN y) { return Fp_mul(x,y,(GEN)E); }
 3690 
 3691 static GEN
 3692 _Fp_sqr(void *E, GEN x) { return Fp_sqr(x,(GEN)E); }
 3693 
 3694 static GEN
 3695 _Fp_one(void *E) { (void) E; return gen_1; }
 3696 
 3697 GEN
 3698 Fp_pow_init(GEN x, GEN n, long k, GEN p)
 3699 {
 3700   return gen_pow_init(x, n, k, (void*)p, &_Fp_sqr, &_Fp_mul);
 3701 }
 3702 
 3703 GEN
 3704 Fp_pow_table(GEN R, GEN n, GEN p)
 3705 {
 3706   return gen_pow_table(R, n, (void*)p, &_Fp_one, &_Fp_mul);
 3707 }
 3708 
 3709 GEN
 3710 Fp_powers(GEN x, long n, GEN p)
 3711 {
 3712   if (lgefint(p) == 3)
 3713     return Flv_to_ZV(Fl_powers(umodiu(x, uel(p, 2)), n, uel(p, 2)));
 3714   return gen_powers(x, n, 1, (void*)p, _Fp_sqr, _Fp_mul, _Fp_one);
 3715 }
 3716 
 3717 GEN
 3718 FpV_prod(GEN V, GEN p)
 3719 {
 3720   return gen_product(V, (void *)p, &_Fp_mul);
 3721 }
 3722 
 3723 static GEN
 3724 _Fp_pow(void *E, GEN x, GEN n) { return Fp_pow(x,n,(GEN)E); }
 3725 
 3726 static GEN
 3727 _Fp_rand(void *E) { return addiu(randomi(subiu((GEN)E,1)),1); }
 3728 
 3729 static GEN Fp_easylog(void *E, GEN a, GEN g, GEN ord);
 3730 
 3731 static const struct bb_group Fp_star={_Fp_mul,_Fp_pow,_Fp_rand,hash_GEN,
 3732                                       equalii,equali1,Fp_easylog};
 3733 
 3734 static GEN
 3735 _Fp_red(void *E, GEN x) { return Fp_red(x, (GEN)E); }
 3736 
 3737 static GEN
 3738 _Fp_add(void *E, GEN x, GEN y) { (void) E; return addii(x,y); }
 3739 
 3740 static GEN
 3741 _Fp_neg(void *E, GEN x) { (void) E; return negi(x); }
 3742 
 3743 static GEN
 3744 _Fp_rmul(void *E, GEN x, GEN y) { (void) E; return mulii(x,y); }
 3745 
 3746 static GEN
 3747 _Fp_inv(void *E, GEN x) { return Fp_inv(x,(GEN)E); }
 3748 
 3749 static int
 3750 _Fp_equal0(GEN x) { return signe(x)==0; }
 3751 
 3752 static GEN
 3753 _Fp_s(void *E, long x) { (void) E; return stoi(x); }
 3754 
 3755 static const struct bb_field Fp_field={_Fp_red,_Fp_add,_Fp_rmul,_Fp_neg,
 3756                                         _Fp_inv,_Fp_equal0,_Fp_s};
 3757 
 3758 const struct bb_field *get_Fp_field(void **E, GEN p)
 3759 {
 3760   *E = (void*)p; return &Fp_field;
 3761 }
 3762 
 3763 /*********************************************************************/
 3764 /**                                                                 **/
 3765 /**               ORDER of INTEGERMOD x  in  (Z/nZ)*                **/
 3766 /**                                                                 **/
 3767 /*********************************************************************/
 3768 ulong
 3769 Fl_order(ulong a, ulong o, ulong p)
 3770 {
 3771   pari_sp av = avma;
 3772   GEN m, P, E;
 3773   long i;
 3774   if (a==1) return 1;
 3775   if (!o) o = p-1;
 3776   m = factoru(o);
 3777   P = gel(m,1);
 3778   E = gel(m,2);
 3779   for (i = lg(P)-1; i; i--)
 3780   {
 3781     ulong j, l = P[i], e = E[i], t = o / upowuu(l,e), y = Fl_powu(a, t, p);
 3782     if (y == 1) o = t;
 3783     else for (j = 1; j < e; j++)
 3784     {
 3785       y = Fl_powu(y, l, p);
 3786       if (y == 1) { o = t *  upowuu(l, j); break; }
 3787     }
 3788   }
 3789   return gc_ulong(av, o);
 3790 }
 3791 
 3792 /*Find the exact order of a assuming a^o==1*/
 3793 GEN
 3794 Fp_order(GEN a, GEN o, GEN p) {
 3795   if (lgefint(p) == 3 && (!o || typ(o) == t_INT))
 3796   {
 3797     ulong pp = p[2], oo = (o && lgefint(o)==3)? uel(o,2): pp-1;
 3798     return utoi( Fl_order(umodiu(a, pp), oo, pp) );
 3799   }
 3800   return gen_order(a, o, (void*)p, &Fp_star);
 3801 }
 3802 GEN
 3803 Fp_factored_order(GEN a, GEN o, GEN p)
 3804 { return gen_factored_order(a, o, (void*)p, &Fp_star); }
 3805 
 3806 /* return order of a mod p^e, e > 0, pe = p^e */
 3807 static GEN
 3808 Zp_order(GEN a, GEN p, long e, GEN pe)
 3809 {
 3810   GEN ap, op;
 3811   if (absequaliu(p, 2))
 3812   {
 3813     if (e == 1) return gen_1;
 3814     if (e == 2) return mod4(a) == 1? gen_1: gen_2;
 3815     if (mod4(a) == 1)
 3816       op = gen_1;
 3817     else {
 3818       op = gen_2;
 3819       a = Fp_sqr(a, pe);
 3820     }
 3821   } else {
 3822     ap = (e == 1)? a: remii(a,p);
 3823     op = Fp_order(ap, subiu(p,1), p);
 3824     if (e == 1) return op;
 3825     a = Fp_pow(a, op, pe); /* 1 mod p */
 3826   }
 3827   if (equali1(a)) return op;
 3828   return mulii(op, powiu(p, e - Z_pval(subiu(a,1), p)));
 3829 }
 3830 
 3831 GEN
 3832 znorder(GEN x, GEN o)
 3833 {
 3834   pari_sp av = avma;
 3835   GEN b, a;
 3836 
 3837   if (typ(x) != t_INTMOD) pari_err_TYPE("znorder [t_INTMOD expected]",x);
 3838   b = gel(x,1); a = gel(x,2);
 3839   if (!equali1(gcdii(a,b))) pari_err_COPRIME("znorder", a,b);
 3840   if (!o)
 3841   {
 3842     GEN fa = Z_factor(b), P = gel(fa,1), E = gel(fa,2);
 3843     long i, l = lg(P);
 3844     o = gen_1;
 3845     for (i = 1; i < l; i++)
 3846     {
 3847       GEN p = gel(P,i);
 3848       long e = itos(gel(E,i));
 3849 
 3850       if (l == 2)
 3851         o = Zp_order(a, p, e, b);
 3852       else {
 3853         GEN pe = powiu(p,e);
 3854         o = lcmii(o, Zp_order(remii(a,pe), p, e, pe));
 3855       }
 3856     }
 3857     return gerepileuptoint(av, o);
 3858   }
 3859   return Fp_order(a, o, b);
 3860 }
 3861 GEN
 3862 order(GEN x) { return znorder(x, NULL); }
 3863 
 3864 /*********************************************************************/
 3865 /**                                                                 **/
 3866 /**               DISCRETE LOGARITHM  in  (Z/nZ)*                   **/
 3867 /**                                                                 **/
 3868 /*********************************************************************/
 3869 static GEN
 3870 Fp_log_halfgcd(ulong bnd, GEN C, GEN g, GEN p)
 3871 {
 3872   pari_sp av = avma;
 3873   GEN h1, h2, F, G;
 3874   if (!Fp_ratlift(g,p,C,shifti(C,-1),&h1,&h2)) return gc_NULL(av);
 3875   if ((F = Z_issmooth_fact(h1, bnd)) && (G = Z_issmooth_fact(h2, bnd)))
 3876   {
 3877     GEN M = cgetg(3, t_MAT);
 3878     gel(M,1) = vecsmall_concat(gel(F, 1),gel(G, 1));
 3879     gel(M,2) = vecsmall_concat(gel(F, 2),zv_neg_inplace(gel(G, 2)));
 3880     return gerepileupto(av, M);
 3881   }
 3882   return gc_NULL(av);
 3883 }
 3884 
 3885 static GEN
 3886 Fp_log_find_rel(GEN b, ulong bnd, GEN C, GEN p, GEN *g, long *e)
 3887 {
 3888   GEN rel;
 3889   do
 3890   {
 3891     (*e)++; *g = Fp_mul(*g, b, p);
 3892     rel = Fp_log_halfgcd(bnd, C, *g, p);
 3893   } while (!rel);
 3894   return rel;
 3895 }
 3896 
 3897 struct Fp_log_rel
 3898 {
 3899   GEN rel;
 3900   ulong prmax;
 3901   long nbrel, nbmax, nbgen;
 3902 };
 3903 
 3904 /* add u^e */
 3905 static void
 3906 addifsmooth1(struct Fp_log_rel *r, GEN z, long u, long e)
 3907 {
 3908   pari_sp av = avma;
 3909   long off = r->prmax+1;
 3910   GEN F = cgetg(3, t_MAT);
 3911   gel(F,1) = vecsmall_append(gel(z,1), off+u);
 3912   gel(F,2) = vecsmall_append(gel(z,2), e);
 3913   gel(r->rel,++r->nbrel) = gerepileupto(av, F);
 3914 }
 3915 
 3916 /* add u^-1 v^-1 */
 3917 static void
 3918 addifsmooth2(struct Fp_log_rel *r, GEN z, long u, long v)
 3919 {
 3920   pari_sp av = avma;
 3921   long off = r->prmax+1;
 3922   GEN P = mkvecsmall2(off+u,off+v), E = mkvecsmall2(-1,-1);
 3923   GEN F = cgetg(3, t_MAT);
 3924   gel(F,1) = vecsmall_concat(gel(z,1), P);
 3925   gel(F,2) = vecsmall_concat(gel(z,2), E);
 3926   gel(r->rel,++r->nbrel) = gerepileupto(av, F);
 3927 }
 3928 
 3929 /*
 3930 Let p=C^2+c
 3931 Solve h = (C+x)*(C+a)-p = 0 [mod l]
 3932 h= -c+x*(C+a)+C*a = 0  [mod l]
 3933 x = (c-C*a)/(C+a) [mod l]
 3934 h = -c+C*(x+a)+a*x
 3935 */
 3936 
 3937 GEN
 3938 Fp_log_sieve_worker(long a, long prmax, GEN C, GEN c, GEN Ci, GEN ci, GEN pi, GEN sz)
 3939 {
 3940   pari_sp ltop = avma;
 3941   long th, n = lg(pi)-1;
 3942   long i, j;
 3943   GEN sieve = zero_zv(a+2)+1;
 3944   GEN L = cgetg(1+a+2, t_VEC);
 3945   pari_sp av = avma;
 3946   long rel = 1;
 3947   GEN z, h;
 3948   h = addis(C,a);
 3949   if ((z = Z_issmooth_fact(h, prmax)))
 3950   {
 3951     gel(L, rel++) = mkvec2(z, mkvecsmall3(1, a, -1));
 3952     av = avma;
 3953   }
 3954   for (i=1; i<=n; i++)
 3955   {
 3956     ulong li = pi[i], s = sz[i], al = a % li;
 3957     ulong u, iv = Fl_invsafe(Fl_add(Ci[i],al,li),li);
 3958     if (!iv) continue;
 3959     u = Fl_mul(Fl_sub(ci[i],Fl_mul(Ci[i],al,li),li), iv ,li);
 3960     for(j = u; j<=a; j+=li)
 3961       sieve[j] += s;
 3962   }
 3963   if (a)
 3964   {
 3965     long e = expi(mulis(C,a));
 3966     th = e - expu(e) - 1;
 3967   } else th = -1;
 3968   for (j=0; j<a; j++)
 3969     if (sieve[j]>=th)
 3970     {
 3971       GEN h = addiu(subii(muliu(C,a+j),c), a*j);
 3972       if ((z = Z_issmooth_fact(h, prmax)))
 3973       {
 3974         gel(L, rel++) = mkvec2(z, mkvecsmall3(2, a, j));
 3975         av = avma;
 3976       } else set_avma(av);
 3977     }
 3978   /* j = a */
 3979   if (sieve[a]>=th)
 3980   {
 3981     GEN h = addiu(subii(muliu(C,2*a),c), a*a);
 3982     if ((z = Z_issmooth_fact(h, prmax)))
 3983       gel(L, rel++) = mkvec2(z, mkvecsmall3(1, a, -2));
 3984   }
 3985   setlg(L, rel);
 3986   return gerepilecopy(ltop, L);
 3987 }
 3988 
 3989 static long
 3990 Fp_log_sieve(struct Fp_log_rel *r, GEN C, GEN c, GEN Ci, GEN ci, GEN pi, GEN sz)
 3991 {
 3992   struct pari_mt pt;
 3993   long i, j, nb = 0;
 3994   GEN worker = snm_closure(is_entry("_Fp_log_sieve_worker"),
 3995                mkvecn(7, utoi(r->prmax), C, c, Ci, ci, pi, sz));
 3996   long running, pending = 0;
 3997   GEN W = zerovec(r->nbgen);
 3998   mt_queue_start_lim(&pt, worker, r->nbgen);
 3999   for (i = 0; (running = (i < r->nbgen)) || pending; i++)
 4000   {
 4001     GEN done;
 4002     long idx;
 4003     mt_queue_submit(&pt, i, running ? mkvec(stoi(i)): NULL);
 4004     done = mt_queue_get(&pt, &idx, &pending);
 4005     if (!done || lg(done)==1) continue;
 4006     gel(W, idx+1) = done;
 4007     nb += lg(done)-1;
 4008     if (DEBUGLEVEL && (i&127)==0)
 4009       err_printf("%ld%% ",100*nb/r->nbmax);
 4010   }
 4011   mt_queue_end(&pt);
 4012   for(j = 1; j <= r->nbgen && r->nbrel < r->nbmax; j++)
 4013   {
 4014     long ll, m;
 4015     GEN L = gel(W,j);
 4016     if (isintzero(L)) continue;
 4017     ll = lg(L);
 4018     for (m=1; m<ll && r->nbrel < r->nbmax ; m++)
 4019     {
 4020       GEN Lm = gel(L,m), h = gel(Lm, 1), v = gel(Lm, 2);
 4021       if (v[1] == 1)
 4022         addifsmooth1(r, h, v[2], v[3]);
 4023       else
 4024         addifsmooth2(r, h, v[2], v[3]);
 4025     }
 4026   }
 4027   return j;
 4028 }
 4029 
 4030 static GEN
 4031 ECP_psi(GEN x, GEN y)
 4032 {
 4033   long prec = realprec(x);
 4034   GEN lx = glog(x, prec), ly = glog(y, prec);
 4035   GEN u = gdiv(lx, ly);
 4036   return gpow(u, gneg(u),prec);
 4037 }
 4038 
 4039 struct computeG
 4040 {
 4041   GEN C;
 4042   long bnd, nbi;
 4043 };
 4044 
 4045 static GEN
 4046 _computeG(void *E, GEN gen)
 4047 {
 4048   struct computeG * d = (struct computeG *) E;
 4049   GEN ps = ECP_psi(gmul(gen,d->C), stoi(d->bnd));
 4050   return gsub(gmul(gsqr(gen),ps),gmul2n(gaddgs(gen,d->nbi),2));
 4051 }
 4052 
 4053 static long
 4054 compute_nbgen(GEN C, long bnd, long nbi)
 4055 {
 4056   struct computeG d;
 4057   d.C = shifti(C, 1);
 4058   d.bnd = bnd;
 4059   d.nbi = nbi;
 4060   return itos(ground(zbrent((void*)&d, _computeG, gen_2, stoi(bnd), DEFAULTPREC)));
 4061 }
 4062 
 4063 static GEN
 4064 _psi(void*E, GEN y)
 4065 {
 4066   GEN lx = (GEN) E;
 4067   long prec = realprec(lx);
 4068   GEN ly = glog(y, prec);
 4069   GEN u = gdiv(lx, ly);
 4070   return gsub(gdiv(y ,ly), gpow(u, u, prec));
 4071 }
 4072 
 4073 static GEN
 4074 opt_param(GEN x, long prec)
 4075 {
 4076   return zbrent((void*)glog(x,prec), _psi, gen_2, x, prec);
 4077 }
 4078 
 4079 static GEN
 4080 check_kernel(long nbg, long N, long prmax, GEN C, GEN M, GEN p, GEN m)
 4081 {
 4082   pari_sp av = avma;
 4083   long lM = lg(M)-1, nbcol = lM;
 4084   long tbs = maxss(1, expu(nbg/expi(m)));
 4085   for (;;)
 4086   {
 4087     GEN K = FpMs_leftkernel_elt_col(M, nbcol, N, m);
 4088     GEN tab;
 4089     long i, f=0;
 4090     long l = lg(K), lm = lgefint(m);
 4091     GEN idx = diviiexact(subiu(p,1),m), g;
 4092     pari_timer ti;
 4093     if (DEBUGLEVEL) timer_start(&ti);
 4094     for(i=1; i<l; i++)
 4095       if (signe(gel(K,i)))
 4096         break;
 4097     g = Fp_pow(utoi(i), idx, p);
 4098     tab = Fp_pow_init(g, p, tbs, p);
 4099     K = FpC_Fp_mul(K, Fp_inv(gel(K,i), m), m);
 4100     for(i=1; i<l; i++)
 4101     {
 4102       GEN k = gel(K,i);
 4103       GEN j = i<=prmax ? utoi(i): addis(C,i-(prmax+1));
 4104       if (signe(k)==0 || !equalii(Fp_pow_table(tab, k, p), Fp_pow(j, idx, p)))
 4105         gel(K,i) = cgetineg(lm);
 4106       else
 4107         f++;
 4108     }
 4109     if (DEBUGLEVEL) timer_printf(&ti,"found %ld/%ld logs", f, nbg);
 4110     if(f > (nbg>>1)) return gerepileupto(av, K);
 4111     for(i=1; i<=nbcol; i++)
 4112     {
 4113       long a = 1+random_Fl(lM);
 4114       swap(gel(M,a),gel(M,i));
 4115     }
 4116     if (4*nbcol>5*nbg) nbcol = nbcol*9/10;
 4117   }
 4118 }
 4119 
 4120 static GEN
 4121 Fp_log_find_ind(GEN a, GEN K, long prmax, GEN C, GEN p, GEN m)
 4122 {
 4123   pari_sp av=avma;
 4124   GEN aa = gen_1;
 4125   long AV = 0;
 4126   for(;;)
 4127   {
 4128     GEN A = Fp_log_find_rel(a, prmax, C, p, &aa, &AV);
 4129     GEN F = gel(A,1), E = gel(A,2);
 4130     GEN Ao = gen_0;
 4131     long i, l = lg(F);
 4132     for(i=1; i<l; i++)
 4133     {
 4134       GEN Ki = gel(K,F[i]);
 4135       if (signe(Ki)<0) break;
 4136       Ao = addii(Ao, mulis(Ki, E[i]));
 4137     }
 4138     if (i==l) return Fp_divu(Ao, AV, m);
 4139     aa = gerepileuptoint(av, aa);
 4140   }
 4141 }
 4142 
 4143 static GEN
 4144 Fp_log_index(GEN a, GEN b, GEN m, GEN p)
 4145 {
 4146   pari_sp av = avma, av2;
 4147   long i, j, nbi, nbr = 0, nbrow, nbg;
 4148   GEN C, c, Ci, ci, pi, pr, sz, l, Ao, Bo, K, d, p_1;
 4149   pari_timer ti;
 4150   struct Fp_log_rel r;
 4151   ulong bnds = itou(roundr_safe(opt_param(sqrti(p),DEFAULTPREC)));
 4152   ulong bnd = 4*bnds;
 4153   if (!bnds || cmpii(sqru(bnds),m)>=0) return NULL;
 4154 
 4155   p_1 = subiu(p,1);
 4156   if (!is_pm1(gcdii(m,diviiexact(p_1,m))))
 4157     m = diviiexact(p_1, Z_ppo(p_1, m));
 4158   pr = primes_upto_zv(bnd);
 4159   nbi = lg(pr)-1;
 4160   C = sqrtremi(p, &c);
 4161   av2 = avma;
 4162   for (i = 1; i <= nbi; ++i)
 4163   {
 4164     ulong lp = pr[i];
 4165     while (lp <= bnd)
 4166     {
 4167       nbr++;
 4168       lp *= pr[i];
 4169     }
 4170   }
 4171   pi = cgetg(nbr+1,t_VECSMALL);
 4172   Ci = cgetg(nbr+1,t_VECSMALL);
 4173   ci = cgetg(nbr+1,t_VECSMALL);
 4174   sz = cgetg(nbr+1,t_VECSMALL);
 4175   for (i = 1, j = 1; i <= nbi; ++i)
 4176   {
 4177     ulong lp = pr[i], sp = expu(2*lp-1);
 4178     while (lp <= bnd)
 4179     {
 4180       pi[j] = lp;
 4181       Ci[j] = umodiu(C, lp);
 4182       ci[j] = umodiu(c, lp);
 4183       sz[j] = sp;
 4184       lp *= pr[i];
 4185       j++;
 4186     }
 4187   }
 4188   r.nbrel = 0;
 4189   r.nbgen = compute_nbgen(C, bnd, nbi);
 4190   r.nbmax = 2*(nbi+r.nbgen);
 4191   r.rel = cgetg(r.nbmax+1,t_VEC);
 4192   r.prmax = pr[nbi];
 4193   if (DEBUGLEVEL)
 4194   {
 4195     err_printf("bnd=%lu Size FB=%ld extra gen=%ld \n", bnd, nbi, r.nbgen);
 4196     timer_start(&ti);
 4197   }
 4198   nbg = Fp_log_sieve(&r, C, c, Ci, ci, pi, sz);
 4199   nbrow = r.prmax + nbg;
 4200   if (DEBUGLEVEL)
 4201   {
 4202     err_printf("\n");
 4203     timer_printf(&ti," %ld relations, %ld generators", r.nbrel, nbi+nbg);
 4204   }
 4205   setlg(r.rel,r.nbrel+1);
 4206   r.rel = gerepilecopy(av2, r.rel);
 4207   K = check_kernel(nbi+nbrow-r.prmax, nbrow, r.prmax, C, r.rel, p, m);
 4208   if (DEBUGLEVEL) timer_start(&ti);
 4209   Ao = Fp_log_find_ind(a, K, r.prmax, C, p, m);
 4210   if (DEBUGLEVEL) timer_printf(&ti," log element");
 4211   Bo = Fp_log_find_ind(b, K, r.prmax, C, p, m);
 4212   if (DEBUGLEVEL) timer_printf(&ti," log generator");
 4213   d = gcdii(Ao,Bo);
 4214   l = Fp_div(diviiexact(Ao, d) ,diviiexact(Bo, d), m);
 4215   if (!equalii(a,Fp_pow(b,l,p))) pari_err_BUG("Fp_log_index");
 4216   return gerepileuptoint(av, l);
 4217 }
 4218 
 4219 static int
 4220 Fp_log_use_index(long e, long p)
 4221 {
 4222   return (e >= 27 && 20*(p+6)<=e*e);
 4223 }
 4224 
 4225 /* Trivial cases a = 1, -1. Return x s.t. g^x = a or [] if no such x exist */
 4226 static GEN
 4227 Fp_easylog(void *E, GEN a, GEN g, GEN ord)
 4228 {
 4229   pari_sp av = avma;
 4230   GEN p = (GEN)E;
 4231   /* assume a reduced mod p, p not necessarily prime */
 4232   if (equali1(a)) return gen_0;
 4233   /* p > 2 */
 4234   if (equalii(subiu(p,1), a))  /* -1 */
 4235   {
 4236     pari_sp av2;
 4237     GEN t;
 4238     ord = get_arith_Z(ord);
 4239     if (mpodd(ord)) { set_avma(av); return cgetg(1, t_VEC); } /* no solution */
 4240     t = shifti(ord,-1); /* only possible solution */
 4241     av2 = avma;
 4242     if (!equalii(Fp_pow(g, t, p), a)) { set_avma(av); return cgetg(1, t_VEC); }
 4243     set_avma(av2); return gerepileuptoint(av, t);
 4244   }
 4245   if (typ(ord)==t_INT && BPSW_psp(p) && Fp_log_use_index(expi(ord),expi(p)))
 4246     return Fp_log_index(a, g, ord, p);
 4247   return gc_NULL(av); /* not easy */
 4248 }
 4249 
 4250 GEN
 4251 Fp_log(GEN a, GEN g, GEN ord, GEN p)
 4252 {
 4253   GEN v = get_arith_ZZM(ord);
 4254   GEN F = gmael(v,2,1);
 4255   long lF = lg(F)-1, lmax;
 4256   if (lF == 0) return equali1(a)? gen_0: cgetg(1, t_VEC);
 4257   lmax = expi(gel(F,lF));
 4258   if (BPSW_psp(p) && Fp_log_use_index(lmax,expi(p)))
 4259     v = mkvec2(gel(v,1),ZM_famat_limit(gel(v,2),int2n(27)));
 4260   return gen_PH_log(a,g,v,(void*)p,&Fp_star);
 4261 }
 4262 
 4263 static ulong
 4264 Fl_log_naive(ulong a, ulong g, ulong ord, ulong p)
 4265 {
 4266   ulong i, h=1;
 4267   for(i=0; i<ord; i++, h = Fl_mul(h, g, p))
 4268     if(a==h) return i;
 4269   return ~0UL;
 4270 }
 4271 
 4272 static ulong
 4273 Fl_log_naive_pre(ulong a, ulong g, ulong ord, ulong p, ulong pi)
 4274 {
 4275   ulong i, h=1;
 4276   for(i=0; i<ord; i++, h = Fl_mul_pre(h, g, p, pi))
 4277     if(a==h) return i;
 4278   return ~0UL;
 4279 }
 4280 
 4281 static ulong
 4282 Fl_log_Fp(ulong a, ulong g, ulong ord, ulong p)
 4283 {
 4284   pari_sp av = avma;
 4285   GEN r = Fp_log(utoi(a),utoi(g),utoi(ord),utoi(p));
 4286   return gc_ulong(av, typ(r)==t_INT ? itou(r): ~0UL);
 4287 }
 4288 
 4289 ulong
 4290 Fl_log_pre(ulong a, ulong g, ulong ord, ulong p, ulong pi)
 4291 {
 4292   if (ord <= 200) return Fl_log_naive_pre(a, g, ord, p, pi);
 4293   return Fl_log_Fp(a, g, ord, p);
 4294 }
 4295 
 4296 ulong
 4297 Fl_log(ulong a, ulong g, ulong ord, ulong p)
 4298 {
 4299   if (ord <= 200)
 4300   return (p&HIGHMASK) ? Fl_log_naive_pre(a, g, ord, p, get_Fl_red(p))
 4301                       : Fl_log_naive(a, g, ord, p);
 4302   return Fl_log_Fp(a, g, ord, p);
 4303 }
 4304 
 4305 /* find x such that h = g^x mod N > 1, N = prod_{i <= l} P[i]^E[i], P[i] prime.
 4306  * PHI[l] = eulerphi(N / P[l]^E[l]).   Destroys P/E */
 4307 static GEN
 4308 znlog_rec(GEN h, GEN g, GEN N, GEN P, GEN E, GEN PHI)
 4309 {
 4310   long l = lg(P) - 1, e = E[l];
 4311   GEN p = gel(P, l), phi = gel(PHI,l), pe = e == 1? p: powiu(p, e);
 4312   GEN a,b, hp,gp, hpe,gpe, ogpe; /* = order(g mod p^e) | p^(e-1)(p-1) */
 4313 
 4314   if (l == 1) {
 4315     hpe = h;
 4316     gpe = g;
 4317   } else {
 4318     hpe = modii(h, pe);
 4319     gpe = modii(g, pe);
 4320   }
 4321   if (e == 1) {
 4322     hp = hpe;
 4323     gp = gpe;
 4324   } else {
 4325     hp = remii(hpe, p);
 4326     gp = remii(gpe, p);
 4327   }
 4328   if (hp == gen_0 || gp == gen_0) return NULL;
 4329   if (absequaliu(p, 2))
 4330   {
 4331     GEN n = int2n(e);
 4332     ogpe = Zp_order(gpe, gen_2, e, n);
 4333     a = Fp_log(hpe, gpe, ogpe, n);
 4334     if (typ(a) != t_INT) return NULL;
 4335   }
 4336   else
 4337   { /* Avoid black box groups: (Z/p^2)^* / (Z/p)^* ~ (Z/pZ, +), where DL
 4338        is trivial */
 4339     /* [order(gp), factor(order(gp))] */
 4340     GEN v = Fp_factored_order(gp, subiu(p,1), p);
 4341     GEN ogp = gel(v,1);
 4342     if (!equali1(Fp_pow(hp, ogp, p))) return NULL;
 4343     a = Fp_log(hp, gp, v, p);
 4344     if (typ(a) != t_INT) return NULL;
 4345     if (e == 1) ogpe = ogp;
 4346     else
 4347     { /* find a s.t. g^a = h (mod p^e), p odd prime, e > 0, (h,p) = 1 */
 4348       /* use p-adic log: O(log p + e) mul*/
 4349       long vpogpe, vpohpe;
 4350 
 4351       hpe = Fp_mul(hpe, Fp_pow(gpe, negi(a), pe), pe);