Code coverage tests

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

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

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

LCOV - code coverage report
Current view: top level - basemath - arith1.c (source / functions) Hit Total Coverage
Test: PARI/GP v2.8.0 lcov report (development 19214-1621e44) Lines: 2617 2847 91.9 %
Date: 2016-07-26 07:10:39 Functions: 229 243 94.2 %
Legend: Lines: hit not hit

          Line data    Source code
       1             : /* Copyright (C) 2000  The PARI group.
       2             : 
       3             : This file is part of the PARI/GP package.
       4             : 
       5             : PARI/GP is free software; you can redistribute it and/or modify it under the
       6             : terms of the GNU General Public License as published by the Free Software
       7             : Foundation. It is distributed in the hope that it will be useful, but WITHOUT
       8             : ANY WARRANTY WHATSOEVER.
       9             : 
      10             : Check the License for details. You should have received a copy of it, along
      11             : with the package; see the file 'COPYING'. If not, write to the Free Software
      12             : Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. */
      13             : 
      14             : /*********************************************************************/
      15             : /**                                                                 **/
      16             : /**                     ARITHMETIC FUNCTIONS                        **/
      17             : /**                         (first part)                            **/
      18             : /**                                                                 **/
      19             : /*********************************************************************/
      20             : #include "pari.h"
      21             : #include "paripriv.h"
      22             : 
      23             : /******************************************************************/
      24             : /*                                                                */
      25             : /*                 GENERATOR of (Z/mZ)*                           */
      26             : /*                                                                */
      27             : /******************************************************************/
      28             : static GEN
      29          42 : remove2(GEN q) { long v = vali(q); return v? shifti(q, -v): q; }
      30             : static ulong
      31       21271 : u_remove2(ulong q) { return q >> vals(q); }
      32             : GEN
      33          42 : odd_prime_divisors(GEN q) { return gel(Z_factor(remove2(q)), 1); }
      34             : static GEN
      35       21271 : u_odd_prime_divisors(ulong q) { return gel(factoru(u_remove2(q)), 1); }
      36             : /* p odd prime, q=(p-1)/2; L0 list of (some) divisors of q = (p-1)/2 or NULL
      37             :  * (all prime divisors of q); return the q/l, l in L0 */
      38             : static GEN
      39          33 : is_gener_expo(GEN p, GEN L0)
      40             : {
      41          33 :   GEN L, q = shifti(p,-1);
      42             :   long i, l;
      43          33 :   if (L0) {
      44          14 :     l = lg(L0);
      45          14 :     L = cgetg(l, t_VEC);
      46             :   } else {
      47          19 :     L0 = L = odd_prime_divisors(q);
      48          19 :     l = lg(L);
      49             :   }
      50          33 :   for (i=1; i<l; i++) gel(L,i) = diviiexact(q, gel(L0,i));
      51          33 :   return L;
      52             : }
      53             : static GEN
      54       21509 : u_is_gener_expo(ulong p, GEN L0)
      55             : {
      56       21509 :   const ulong q = p >> 1;
      57             :   long i, l;
      58             :   GEN L;
      59       21509 :   if (L0) {
      60        2219 :     l = lg(L0);
      61        2219 :     L = cgetg(l, t_VECSMALL);
      62             :   } else {
      63       19290 :     L0 = L = u_odd_prime_divisors(q);
      64       19290 :     l = lg(L);
      65             :   }
      66       21509 :   for (i=1; i<l; i++) L[i] = q / uel(L0,i);
      67       21509 :   return L;
      68             : }
      69             : 
      70             : int
      71       53336 : is_gener_Fl(ulong x, ulong p, ulong p_1, GEN L)
      72             : {
      73             :   long i;
      74       53336 :   if (krouu(x, p) >= 0) return 0;
      75       68417 :   for (i=lg(L)-1; i; i--)
      76             :   {
      77       45468 :     ulong t = Fl_powu(x, uel(L,i), p);
      78       45468 :     if (t == p_1 || t == 1) return 0;
      79             :   }
      80       22949 :   return 1;
      81             : }
      82             : /* assume p prime */
      83             : ulong
      84       51860 : pgener_Fl_local(ulong p, GEN L0)
      85             : {
      86       51860 :   const pari_sp av = avma;
      87       51860 :   const ulong p_1 = p-1;
      88             :   long x;
      89             :   GEN L;
      90       51860 :   if (p <= 19) switch(p)
      91             :   { /* quick trivial cases */
      92          42 :     case 2:  return 1;
      93             :     case 7:
      94        3845 :     case 17: return 3;
      95       26488 :     default: return 2;
      96             :   }
      97       21485 :   L = u_is_gener_expo(p,L0);
      98       21485 :   for (x=2;;x++) { if (is_gener_Fl(x,p,p_1,L)) break; }
      99       21485 :   avma = av; return x;
     100             : }
     101             : ulong
     102       43594 : pgener_Fl(ulong p) { return pgener_Fl_local(p, NULL); }
     103             : 
     104             : /* L[i] = set of (p-1)/2l, l ODD prime divisor of p-1 (l=2 can be included,
     105             :  * but wasteful) */
     106             : int
     107         141 : is_gener_Fp(GEN x, GEN p, GEN p_1, GEN L)
     108             : {
     109         141 :   long i, t = lgefint(x)==3? kroui(x[2], p): kronecker(x, p);
     110         141 :   if (t >= 0) return 0;
     111         381 :   for (i = lg(L)-1; i; i--)
     112             :   {
     113         255 :     GEN t = Fp_pow(x, gel(L,i), p);
     114         255 :     if (equalii(t, p_1) || equali1(t)) return 0;
     115             :   }
     116         126 :   return 1;
     117             : }
     118             : 
     119             : /* assume p prime, return a generator of all L[i]-Sylows in F_p^*. */
     120             : GEN
     121        7574 : pgener_Fp_local(GEN p, GEN L0)
     122             : {
     123        7574 :   pari_sp av0 = avma;
     124             :   GEN x, p_1, L;
     125        7574 :   if (lgefint(p) == 3)
     126             :   {
     127             :     ulong z;
     128        7545 :     if (p[2] == 2) return gen_1;
     129        4780 :     if (L0) L0 = ZV_to_nv(L0);
     130        4780 :     z = pgener_Fl_local(uel(p,2), L0);
     131        4780 :     avma = av0; return utoipos(z);
     132             :   }
     133          29 :   p_1 = subis(p,1); L = is_gener_expo(p, L0);
     134          29 :   x = utoipos(2);
     135          32 :   for (;; x[2]++) { if (is_gener_Fp(x, p, p_1, L)) break; }
     136          29 :   avma = av0; return utoipos(uel(x,2));
     137             : }
     138             : 
     139             : GEN
     140        7042 : pgener_Fp(GEN p) { return pgener_Fp_local(p, NULL); }
     141             : 
     142             : ulong
     143        3797 : pgener_Zl(ulong p)
     144             : {
     145        3797 :   if (p == 2) pari_err_DOMAIN("pgener_Zl","p","=",gen_2,gen_2);
     146             :   /* only p < 2^32 such that znprimroot(p) != znprimroot(p^2) */
     147        3797 :   if (p == 40487) return 10;
     148             : #ifndef LONG_IS_64BIT
     149         539 :   return pgener_Fl(p);
     150             : #else
     151        3258 :   if (p < (1UL<<32)) return pgener_Fl(p);
     152             :   else
     153             :   {
     154          24 :     const pari_sp av = avma;
     155          24 :     const ulong p_1 = p-1;
     156             :     long x ;
     157          24 :     GEN p2 = sqru(p), L = u_is_gener_expo(p, NULL);
     158          96 :     for (x=2;;x++)
     159          96 :       if (is_gener_Fl(x,p,p_1,L) && !is_pm1(Fp_powu(utoipos(x),p_1,p2))) break;
     160          72 :     avma = av; return x;
     161             :   }
     162             : #endif
     163             : }
     164             : 
     165             : /* p prime. Return a primitive root modulo p^e, e > 1 */
     166             : GEN
     167        3801 : pgener_Zp(GEN p)
     168             : {
     169        3801 :   if (lgefint(p) == 3) return utoipos(pgener_Zl(p[2]));
     170             :   else
     171             :   {
     172           4 :     const pari_sp av = avma;
     173           4 :     GEN p_1 = subis(p,1), p2 = sqri(p), L = is_gener_expo(p,NULL);
     174           4 :     GEN x = utoipos(2);
     175          12 :     for (;; x[2]++)
     176          16 :       if (is_gener_Fp(x,p,p_1,L) && !equali1(Fp_pow(x,p_1,p2))) break;
     177          12 :     avma = av; return utoipos(uel(x,2));
     178             :   }
     179             : }
     180             : 
     181             : static GEN
     182         287 : gener_Zp(GEN q, GEN F)
     183             : {
     184         287 :   GEN p = NULL;
     185         287 :   long e = 0;
     186         287 :   if (F)
     187             :   {
     188          14 :     GEN P = gel(F,1), E = gel(F,2);
     189          14 :     long i, l = lg(P);
     190          42 :     for (i = 1; i < l; i++)
     191             :     {
     192          28 :       p = gel(P,i);
     193          28 :       if (equaliu(p, 2)) continue;
     194          14 :       if (i < l-1) pari_err_DOMAIN("znprimroot", "argument","=",F,F);
     195          14 :       e = itos(gel(E,i));
     196             :     }
     197          14 :     if (!p) pari_err_DOMAIN("znprimroot", "argument","=",F,F);
     198             :   }
     199             :   else
     200         273 :     e = Z_isanypower(q, &p);
     201         287 :   return e > 1? pgener_Zp(p): pgener_Fp(q);
     202             : }
     203             : 
     204             : GEN
     205         357 : znprimroot(GEN N)
     206             : {
     207         357 :   pari_sp av = avma;
     208             :   GEN x, n, F;
     209             : 
     210         357 :   if ((F = check_arith_non0(N,"znprimroot")))
     211             :   {
     212          14 :     F = clean_Z_factor(F);
     213          14 :     N = typ(N) == t_VEC? gel(N,1): factorback(F);
     214             :   }
     215         350 :   if (signe(N) < 0) N = absi(N);
     216         350 :   if (cmpiu(N, 4) <= 0) { avma = av; return mkintmodu(N[2]-1,N[2]); }
     217         301 :   switch(mod4(N))
     218             :   {
     219             :     case 0: /* N = 0 mod 4 */
     220          14 :       pari_err_DOMAIN("znprimroot", "argument","=",N,N);
     221           0 :       x = NULL; break;
     222             :     case 2: /* N = 2 mod 4 */
     223          21 :       n = shifti(N,-1); /* becomes odd */
     224          21 :       x = gener_Zp(n,F); if (!mod2(x)) x = addii(x,n);
     225          21 :       break;
     226             :     default: /* N odd */
     227         266 :       x = gener_Zp(N,F);
     228         266 :       break;
     229             :   }
     230         287 :   return gerepilecopy(av, mkintmod(x, N));
     231             : }
     232             : 
     233             : /* n | (p-1), returns a primitive n-th root of 1 in F_p^* */
     234             : GEN
     235           0 : rootsof1_Fp(GEN n, GEN p)
     236             : {
     237           0 :   pari_sp av = avma;
     238           0 :   GEN L = odd_prime_divisors(n); /* 2 implicit in pgener_Fp_local */
     239           0 :   GEN z = pgener_Fp_local(p, L);
     240           0 :   z = Fp_pow(z, diviiexact(subis(p,1), n), p); /* prim. n-th root of 1 */
     241           0 :   return gerepileuptoint(av, z);
     242             : }
     243             : 
     244             : GEN
     245         301 : rootsof1u_Fp(ulong n, GEN p)
     246             : {
     247         301 :   pari_sp av = avma;
     248         301 :   GEN z, L = u_odd_prime_divisors(n); /* 2 implicit in pgener_Fp_local */
     249         301 :   z = pgener_Fp_local(p, Flv_to_ZV(L));
     250         301 :   z = Fp_pow(z, diviuexact(subis(p,1), n), p); /* prim. n-th root of 1 */
     251         301 :   return gerepileuptoint(av, z);
     252             : }
     253             : 
     254             : ulong
     255        1680 : rootsof1_Fl(ulong n, ulong p)
     256             : {
     257        1680 :   pari_sp av = avma;
     258        1680 :   GEN L = u_odd_prime_divisors(n); /* 2 implicit in pgener_Fl_local */
     259        1680 :   ulong z = pgener_Fl_local(p, L);
     260        1680 :   z = Fl_powu(z, (p-1) / n, p); /* prim. n-th root of 1 */
     261        1680 :   avma = av; return z;
     262             : }
     263             : 
     264             : /*********************************************************************/
     265             : /**                                                                 **/
     266             : /**                     INVERSE TOTIENT FUNCTION                    **/
     267             : /**                                                                 **/
     268             : /*********************************************************************/
     269             : /* N t_INT, L a ZV containing all prime divisors of N, and possibly other
     270             :  * primes. Return factor(N) */
     271             : GEN
     272      350651 : Z_factor_listP(GEN N, GEN L)
     273             : {
     274      350651 :   long i, k, l = lg(L);
     275      350651 :   GEN P = cgetg(l, t_COL), E = cgetg(l, t_COL);
     276     1346688 :   for (i = k = 1; i < l; i++)
     277             :   {
     278      996037 :     GEN p = gel(L,i);
     279      996037 :     long v = Z_pvalrem(N, p, &N);
     280      996037 :     if (v)
     281             :     {
     282      792176 :       gel(P,k) = p;
     283      792176 :       gel(E,k) = utoipos(v);
     284      792176 :       k++;
     285             :     }
     286             :   }
     287      350651 :   setlg(P, k);
     288      350651 :   setlg(E, k); return mkmat2(P,E);
     289             : }
     290             : 
     291             : /* look for x such that phi(x) = n, p | x => p > m (if m = NULL: no condition).
     292             :  * L is a list of primes containing all prime divisors of n. */
     293             : static long
     294      621565 : istotient_i(GEN n, GEN m, GEN L, GEN *px)
     295             : {
     296      621565 :   pari_sp av = avma, av2;
     297             :   GEN k, D;
     298             :   long i, v;
     299      621565 :   if (m && mod2(n))
     300             :   {
     301      270914 :     if (!equali1(n)) return 0;
     302       69986 :     if (px) *px = gen_1;
     303       69986 :     return 1;
     304             :   }
     305      350651 :   D = divisors(Z_factor_listP(shifti(n, -1), L));
     306             :   /* loop through primes p > m, d = p-1 | n */
     307      350651 :   av2 = avma;
     308      350651 :   if (!m)
     309             :   { /* special case p = 2, d = 1 */
     310       69986 :     k = n;
     311       69986 :     for (v = 1;; v++) {
     312       69986 :       if (istotient_i(k, gen_2, L, px)) {
     313       69986 :         if (px) *px = shifti(*px, v);
     314       69986 :         return 1;
     315             :       }
     316           0 :       if (mod2(k)) break;
     317           0 :       k = shifti(k,-1);
     318           0 :     }
     319           0 :     avma = av2;
     320             :   }
     321     1099462 :   for (i = 1; i < lg(D); ++i)
     322             :   {
     323     1001588 :     GEN p, d = shifti(gel(D, i), 1); /* even divisors of n */
     324     1001588 :     if (m && cmpii(d, m) < 0) continue;
     325      677782 :     p = addiu(d, 1);
     326      677782 :     if (!isprime(p)) continue;
     327      442064 :     k = diviiexact(n, d);
     328      481593 :     for (v = 1;; v++) {
     329             :       GEN r;
     330      481593 :       if (istotient_i(k, p, L, px)) {
     331      182791 :         if (px) *px = mulii(*px, powiu(p, v));
     332      182791 :         return 1;
     333             :       }
     334      298802 :       k = dvmdii(k, p, &r);
     335      298802 :       if (r != gen_0) break;
     336       39529 :     }
     337      259273 :     avma = av2;
     338             :   }
     339       97874 :   avma = av; return 0;
     340             : }
     341             : 
     342             : /* find x such that phi(x) = n */
     343             : long
     344       70000 : istotient(GEN n, GEN *px)
     345             : {
     346       70000 :   pari_sp av = avma;
     347       70000 :   if (typ(n) != t_INT) pari_err_TYPE("istotient", n);
     348       70000 :   if (signe(n) < 1) return 0;
     349       70000 :   if (mod2(n))
     350             :   {
     351          14 :     if (!equali1(n)) return 0;
     352          14 :     if (px) *px = gen_1;
     353          14 :     return 1;
     354             :   }
     355       69986 :   if (istotient_i(n, NULL, gel(Z_factor(n), 1), px))
     356             :   {
     357       69986 :     if (!px) avma = av;
     358             :     else
     359       69986 :       *px = gerepileuptoint(av, *px);
     360       69986 :     return 1;
     361             :   }
     362           0 :   avma = av; return 0;
     363             : }
     364             : 
     365             : /*********************************************************************/
     366             : /**                                                                 **/
     367             : /**                     INTEGRAL LOGARITHM                          **/
     368             : /**                                                                 **/
     369             : /*********************************************************************/
     370             : 
     371             : /* y > 1 and B > 0 integers. Return e such that y^(e-1) <= B < y^e, i.e
     372             :  * e = 1 + floor(log_y B). Set *ptq = y^e if non-NULL */
     373             : long
     374      309519 : logint(GEN B, GEN y, GEN *ptq)
     375             : {
     376      309519 :   pari_sp av = avma, av2;
     377             :   long eB, ey, e, i, fl;
     378      309519 :   GEN q,pow2, r = y;
     379             : 
     380      309519 :   if (typ(B) != t_INT) B = ceil_safe(B);
     381      309519 :   eB = expi(B); /* 2^eB <= B < 2^(eB + 1) */
     382      309519 :   ey = expi(y); /* result e satisfies e > eB / (ey+1) */
     383      309519 :   if (eB <= 13 * ey) /* e small, be naive */
     384             :   {
     385      783202 :     for (e=1;; e++)
     386             :     { /* here, r = y^e */
     387      783202 :       fl = cmpii(r, B);
     388      783202 :       if (fl > 0) goto END;
     389      496567 :       r = mulii(r,y);
     390      496567 :     }
     391             :   }
     392             :   /* e > 13 ey / (ey + 1) >= 6.5 */
     393             : 
     394             :   /* binary splitting: compute bits of e one by one */
     395             :   /* compute pow2[i] = y^(2^i) [i < crude upper bound for log_2 log_y(B)] */
     396       22884 :   pow2 = new_chunk((long)log2(eB)+2);
     397       22884 :   gel(pow2,0) = y;
     398       22884 :   for (i=0,q=r;; )
     399             :   { /* r = y^2^i */
     400      124903 :     fl = cmpii(r,B);
     401      124903 :     if (!fl) { e = 1L<<i; e++; r = mulii(r,y); goto END; }
     402      124903 :     if (fl > 0) break;
     403      102019 :     q = r; r = sqri(q);
     404      102019 :     i++; gel(pow2,i) = r;
     405      102019 :   }
     406             : 
     407       22884 :   av2 = avma;
     408       22884 :   for (i--, e=1L<<i;;)
     409             :   { /* y^e = q < B <= r = q * y^(2^i) */
     410      102005 :     if (!fl) break; /* B = r */
     411             :     /* q < B < r */
     412      101998 :     if (--i < 0) { if (fl > 0) e++; break; }
     413       79121 :     r = mulii(q, gel(pow2,i));
     414       79121 :     fl = cmpii(r, B);
     415       79121 :     if (fl <= 0) { e += (1L<<i); q = r = gerepileuptoint(av2, r); }
     416       79121 :   }
     417       22884 :   if (fl <= 0) { e++; r = mulii(r,y); }
     418             : END:
     419      309519 :   if (ptq) *ptq = gerepileuptoint(av, r); else avma = av;
     420      309519 :   return e;
     421             : }
     422             : 
     423             : long
     424          42 : logint0(GEN B, GEN y, GEN *ptq)
     425             : {
     426             :   long e;
     427          42 :   if (typ(B) != t_INT) pari_err_TYPE("logint",B);
     428          42 :   if (signe(B)<=0)
     429           0 :     pari_err_DOMAIN("logint", "x" ,"<=", gen_0, B);
     430          42 :   if (typ(y) != t_INT) pari_err_TYPE("logint",y);
     431          42 :   if (signe(y)<=0 || equali1(y))
     432           0 :     pari_err_DOMAIN("logint", "b" ,"<=", gen_1, y);
     433          42 :   if (equaliu(y, 2))
     434             :   {
     435          21 :     e = expi(B);
     436          21 :     if (ptq) *ptq = int2n(e);
     437          21 :     return e;
     438             :   }
     439          21 :   e = logint(B,y,ptq)-1;
     440          21 :   if (ptq)
     441             :   {
     442           0 :     pari_sp av = avma;
     443           0 :     *ptq = gerepileuptoint(av, diviiexact(*ptq, y));
     444             :   }
     445          21 :   return e;
     446             : }
     447             : 
     448             : /*********************************************************************/
     449             : /**                                                                 **/
     450             : /**                     INTEGRAL SQUARE ROOT                        **/
     451             : /**                                                                 **/
     452             : /*********************************************************************/
     453             : GEN
     454       38669 : sqrtint(GEN a)
     455             : {
     456       38669 :   if (typ(a) != t_INT) pari_err_TYPE("sqrtint",a);
     457       38669 :   switch (signe(a))
     458             :   {
     459       38655 :     case 1: return sqrti(a);
     460           7 :     case 0: return gen_0;
     461           7 :     default: pari_err_DOMAIN("sqrtint", "argument", "<", gen_0,a);
     462             :   }
     463           0 :   return NULL; /* not reached */
     464             : }
     465             : 
     466             : /*********************************************************************/
     467             : /**                                                                 **/
     468             : /**                      PERFECT SQUARE                             **/
     469             : /**                                                                 **/
     470             : /*********************************************************************/
     471             : static int
     472     8773769 : carremod(ulong A)
     473             : {
     474     8773769 :   const int carresmod64[]={
     475             :     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,
     476             :     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};
     477     8773769 :   const int carresmod63[]={
     478             :     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,
     479             :     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};
     480     8773769 :   const int carresmod65[]={
     481             :     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,
     482             :     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};
     483     8773769 :   const int carresmod11[]={1,1,0,1,1,1,0,0,0,1, 0};
     484    17547538 :   return (carresmod64[A & 0x3fUL]
     485     3229936 :     && carresmod63[A % 63UL]
     486     1939294 :     && carresmod65[A % 65UL]
     487    10361511 :     && carresmod11[A % 11UL]);
     488             : }
     489             : 
     490             : /* emulate Z_issquareall on single-word integers */
     491             : long
     492     8607641 : uissquareall(ulong A, ulong *sqrtA)
     493             : {
     494     8607641 :   if (!A) { *sqrtA = 0; return 1; }
     495     8607641 :   if (carremod(A))
     496             :   {
     497     1497173 :     ulong a = usqrt(A);
     498     1497154 :     if (a * a == A) { *sqrtA = a; return 1; }
     499             :   }
     500     7200427 :   return 0;
     501             : }
     502             : long
     503      115450 : uissquare(ulong A)
     504             : {
     505      115450 :   if (!A) return 1;
     506      115450 :   if (carremod(A))
     507             :   {
     508        2335 :     ulong a = usqrt(A);
     509        2335 :     if (a * a == A) return 1;
     510             :   }
     511      113139 :   return 0;
     512             : }
     513             : 
     514             : long
     515     4605926 : Z_issquareall(GEN x, GEN *pt)
     516             : {
     517             :   pari_sp av;
     518             :   GEN y, r;
     519             : 
     520     4605926 :   switch(signe(x))
     521             :   {
     522     1856209 :     case -1: return 0;
     523         686 :     case 0: if (pt) *pt=gen_0; return 1;
     524             :   }
     525     2749031 :   if (lgefint(x) == 3)
     526             :   {
     527     2698513 :     ulong u = uel(x,2), a;
     528     2698513 :     if (!pt) return uissquare(u);
     529     2583063 :     if (!uissquareall(u, &a)) return 0;
     530     1376179 :     *pt = utoipos(a); return 1;
     531             :   }
     532       50518 :   if (!carremod(umodiu(x, 64*63*65*11))) return 0;
     533        6263 :   av = avma; y = sqrtremi(x, &r);
     534        6263 :   if (r != gen_0) { avma = av; return 0; }
     535        5676 :   if (pt) { *pt = y; avma = (pari_sp)y; } else avma = av;
     536        5676 :   return 1;
     537             : }
     538             : 
     539             : /* a t_INT, p prime */
     540             : long
     541           0 : Zp_issquare(GEN a, GEN p)
     542             : {
     543             :   long v;
     544             :   GEN ap;
     545             : 
     546           0 :   if (!signe(a) || gequal1(a)) return 1;
     547           0 :   v = Z_pvalrem(a, p, &ap);
     548           0 :   if (v&1) return 0;
     549           0 :   return equaliu(p, 2)? umodiu(ap, 8) == 1
     550           0 :                       : kronecker(ap,p) == 1;
     551             : }
     552             : 
     553             : static long
     554        1603 : polissquareall(GEN x, GEN *pt)
     555             : {
     556             :   pari_sp av;
     557             :   long v;
     558             :   GEN y, a, b, p;
     559             : 
     560        1603 :   if (!signe(x))
     561             :   {
     562           7 :     if (pt) *pt = gcopy(x);
     563           7 :     return 1;
     564             :   }
     565        1596 :   if (odd(degpol(x))) return 0; /* odd degree */
     566        1057 :   av = avma;
     567        1057 :   v = RgX_valrem(x, &x);
     568        1057 :   if (v & 1) { avma = av; return 0; }
     569        1050 :   a = gel(x,2); /* test constant coeff */
     570        1050 :   if (!pt)
     571          63 :   { if (!issquare(a)) { avma = av; return 0; } }
     572             :   else
     573         987 :   { if (!issquareall(a,&b)) { avma = av; return 0; } }
     574        1050 :   if (!degpol(x)) { /* constant polynomial */
     575          77 :     if (!pt) { avma = av; return 1; }
     576          35 :     y = scalarpol(b, varn(x)); goto END;
     577             :   }
     578         973 :   p = characteristic(x);
     579         973 :   if (signe(p) && !mod2(p))
     580             :   {
     581             :     long i, lx;
     582          35 :     if (!equaliu(p,2)) pari_err_IMPL("issquare for even characteristic != 2");
     583          28 :     x = gmul(x, mkintmod(gen_1, gen_2));
     584          28 :     lx = lg(x);
     585          28 :     if ((lx-3) & 1) { avma = av; return 0; }
     586          49 :     for (i = 3; i < lx; i+=2)
     587          28 :       if (!gequal0(gel(x,i))) { avma = av; return 0; }
     588          21 :     if (pt) {
     589          14 :       y = cgetg((lx+3) / 2, t_POL);
     590          49 :       for (i = 2; i < lx; i+=2)
     591          35 :         if (!issquareall(gel(x,i), &gel(y, (i+2)>>1))) { avma = av; return 0; }
     592          14 :       y[1] = evalsigne(1) | evalvarn(varn(x));
     593          14 :       goto END;
     594             :     } else {
     595          21 :       for (i = 2; i < lx; i+=2)
     596          14 :         if (!issquare(gel(x,i))) { avma = av; return 0; }
     597           7 :       avma = av; return 1;
     598             :     }
     599             :   }
     600             :   else
     601             :   {
     602         938 :     long m = 1;
     603         938 :     x = RgX_Rg_div(x,a);
     604             :     /* a(x^m) = B^2 => B = b(x^m) provided a(0) != 0 */
     605         938 :     if (!signe(p)) x = RgX_deflate_max(x,&m);
     606         938 :     y = ser2rfrac_i(gsqrt(RgX_to_ser(x,lg(x)-1),0));
     607        1190 :     if (!RgX_equal(RgX_sqr(y), x)) { avma = av; return 0; }
     608         686 :     if (!pt) { avma = av; return 1; }
     609         686 :     if (!gequal1(a)) y = gmul(b, y);
     610         686 :     if (m != 1) y = RgX_inflate(y,m);
     611             :   }
     612             : END:
     613         735 :   if (v) y = RgX_shift_shallow(y, v>>1);
     614         735 :   *pt = gerepilecopy(av, y); return 1;
     615             : }
     616             : 
     617             : /* b unit mod p */
     618             : static int
     619         294 : Up_ispower(GEN b, GEN K, GEN p, long d, GEN *pt)
     620             : {
     621         294 :   if (d == 1)
     622             :   { /* mod p: faster */
     623         210 :     if (!Fp_ispower(b, K, p)) return 0;
     624         210 :     if (pt) *pt = Fp_sqrtn(b, K, p, NULL);
     625             :   }
     626             :   else
     627             :   { /* mod p^{2 +} */
     628          84 :     if (!ispower(cvtop(b, p, d), K, pt)) return 0;
     629          63 :     if (pt) *pt = gtrunc(*pt);
     630             :   }
     631         273 :   return 1;
     632             : }
     633             : 
     634             : /* We're studying whether a mod (q*p^e) is a K-th power, (q,p) = 1.
     635             :  * Decide mod p^e, then reduce a mod q unless q = NULL. */
     636             : static int
     637         434 : handle_pe(GEN *pa, GEN q, GEN L, GEN K, GEN p, long e)
     638             : {
     639             :   GEN t, A;
     640         434 :   long v = Z_pvalrem(*pa, p, &A), d = e - v;
     641         434 :   if (d <= 0) t = gen_0;
     642             :   else
     643             :   {
     644             :     ulong r;
     645         378 :     v = udivui_rem(v, K, &r);
     646         378 :     if (r || !Up_ispower(A, K, p, d, L? &t: NULL)) return 0;
     647         273 :     if (L && v) t = mulii(t, powiu(p, v));
     648             :   }
     649         329 :   if (q) *pa = modii(*pa, q);
     650         329 :   if (L) vectrunc_append(L, mkintmod(t, powiu(p, e)));
     651         329 :   return 1;
     652             : }
     653             : long
     654         322 : Zn_ispower(GEN a, GEN q, GEN K, GEN *pt)
     655             : {
     656             :   GEN L, N;
     657             :   pari_sp av;
     658             :   long e, i, l;
     659             :   ulong pp;
     660             :   forprime_t S;
     661             : 
     662         322 :   if (!signe(a))
     663             :   {
     664           7 :     if (pt) {
     665           7 :       GEN t = cgetg(3, t_INTMOD);
     666           7 :       gel(t,1) = icopy(q); gel(t,2) = gen_0; *pt = t;
     667             :     }
     668           7 :     return 1;
     669             :   }
     670             :   /* a != 0 */
     671         315 :   av = avma;
     672             : 
     673         315 :   if (typ(q) != t_INT) /* integer factorization */
     674             :   {
     675           0 :     GEN P = gel(q,1), E = gel(q,2);
     676           0 :     l = lg(P);
     677           0 :     L = pt? vectrunc_init(l): NULL;
     678           0 :     for (i = 1; i < l; i++)
     679             :     {
     680           0 :       GEN p = gel(P,i);
     681           0 :       long e = itos(gel(E,i));
     682           0 :       if (!handle_pe(&a, NULL, L, K, p, e)) { avma = av; return 0; }
     683             :     }
     684           0 :     goto END;
     685             :   }
     686         315 :   if (!mod2(K)
     687         189 :       && kronecker(a, shifti(q,-vali(q))) == -1) { avma = av; return 0; }
     688         308 :   L = pt? vectrunc_init(expi(q)+1): NULL;
     689         308 :   u_forprime_init(&S, 2, tridiv_bound(q));
     690         308 :   while ((pp = u_forprime_next(&S)))
     691             :   {
     692             :     int stop;
     693      883204 :     e = Z_lvalrem_stop(&q, pp, &stop);
     694      883204 :     if (!e) continue;
     695         210 :     if (!handle_pe(&a, q, L, K, utoipos(pp), e)) { avma = av; return 0; }
     696         168 :     if (stop)
     697             :     {
     698         133 :       if (!is_pm1(q) && !handle_pe(&a, q, L, K, q, 1)) { avma = av; return 0; }
     699         133 :       goto END;
     700             :     }
     701             :   }
     702         154 :   l = lg(primetab);
     703         154 :   for (i = 1; i < l; i++)
     704             :   {
     705           0 :     GEN p = gel(primetab,i);
     706           0 :     e = Z_pvalrem(q, p, &q);
     707           0 :     if (!e) continue;
     708           0 :     if (!handle_pe(&a, q, L, K, p, e)) { avma = av; return 0; }
     709           0 :     if (is_pm1(q)) goto END;
     710             :   }
     711         154 :   N = gcdii(a,q);
     712         154 :   if (!is_pm1(N))
     713             :   {
     714         112 :     if (ifac_isprime(N))
     715             :     {
     716          70 :       e = Z_pvalrem(q, N, &q);
     717          70 :       if (!handle_pe(&a, q, L, K, N, e)) { avma = av; return 0; }
     718             :     }
     719             :     else
     720             :     {
     721          42 :       GEN part = ifac_start(N, 0);
     722             :       for(;;)
     723             :       {
     724             :         long e;
     725             :         GEN p;
     726          84 :         if (!ifac_next(&part, &p, &e)) break;
     727          42 :         e = Z_pvalrem(q, p, &q);
     728          42 :         if (!handle_pe(&a, q, L, K, p, e)) { avma = av; return 0; }
     729          42 :       }
     730             :     }
     731             :   }
     732          84 :   if (!is_pm1(q))
     733             :   {
     734          84 :     if (ifac_isprime(q))
     735             :     {
     736          28 :       if (!handle_pe(&a, q, L, K, q, 1)) { avma = av; return 0; }
     737             :     }
     738             :     else
     739             :     {
     740          56 :       GEN part = ifac_start(q, 0);
     741             :       for(;;)
     742             :       {
     743             :         long e;
     744             :         GEN p;
     745         140 :         if (!ifac_next(&part, &p, &e)) break;
     746          98 :         if (!handle_pe(&a, q, L, K, p, e)) { avma = av; return 0; }
     747          84 :       }
     748             :     }
     749             :   }
     750             : END:
     751         203 :   if (pt) *pt = gerepileupto(av, chinese1_coprime_Z(L));
     752         203 :   return 1;
     753             : }
     754             : 
     755             : static long
     756         126 : polmodispower(GEN x, GEN K, GEN *pt)
     757             : {
     758         126 :   pari_sp av = avma;
     759         126 :   GEN p = NULL, T = NULL;
     760         126 :   if (Rg_is_FpXQ(x, &T,&p) && p)
     761             :   {
     762         112 :     x = liftall_shallow(x);
     763         112 :     if (!Fq_ispower(x, K, T, p)) { avma = av; return 0; }
     764         105 :     if (!pt) { avma = av; return 1; }
     765          98 :     x = Fq_sqrtn(x, K, T,p, NULL);
     766          98 :     if (typ(x) == t_INT)
     767           7 :       x = Fp_to_mod(x,p);
     768             :     else
     769          91 :       x = mkpolmod(FpX_to_mod(x,p), FpX_to_mod(T,p));
     770          98 :     *pt = gerepilecopy(av, x); return 1;
     771             :   }
     772          14 :   pari_err_IMPL("ispower for general t_POLMOD");
     773           0 :   return 0;
     774             : }
     775             : 
     776             : long
     777      146174 : issquareall(GEN x, GEN *pt)
     778             : {
     779      146174 :   long tx = typ(x);
     780             :   GEN F;
     781             :   pari_sp av;
     782             : 
     783      146174 :   if (!pt) return issquare(x);
     784        3339 :   switch(tx)
     785             :   {
     786        1442 :     case t_INT: return Z_issquareall(x, pt);
     787         238 :     case t_FRAC: av = avma;
     788         238 :       F = cgetg(3, t_FRAC);
     789         238 :       if (   !Z_issquareall(gel(x,1), &gel(F,1))
     790         105 :           || !Z_issquareall(gel(x,2), &gel(F,2))) { avma = av; return 0; }
     791         105 :       *pt = F; return 1;
     792             : 
     793             :     case t_POLMOD:
     794          21 :       return polmodispower(x, gen_2, pt);
     795        1526 :     case t_POL: return polissquareall(x,pt);
     796           7 :     case t_RFRAC: av = avma;
     797           7 :       F = cgetg(3, t_RFRAC);
     798           7 :       if (   !issquareall(gel(x,1), &gel(F,1))
     799           7 :           || !polissquareall(gel(x,2), &gel(F,2))) { avma = av; return 0; }
     800           7 :       *pt = F; return 1;
     801             : 
     802             :     case t_REAL: case t_COMPLEX: case t_PADIC: case t_SER:
     803           7 :       if (!issquare(x)) return 0;
     804           7 :       *pt = gsqrt(x, DEFAULTPREC); return 1;
     805             : 
     806             :     case t_INTMOD:
     807          63 :       return Zn_ispower(gel(x,2), gel(x,1), gen_2, pt);
     808             : 
     809          35 :     case t_FFELT: return FF_issquareall(x, pt);
     810             : 
     811             :   }
     812           0 :   pari_err_TYPE("issquareall",x);
     813           0 :   return 0; /* not reached */
     814             : }
     815             : 
     816             : long
     817      143024 : issquare(GEN x)
     818             : {
     819             :   pari_sp av;
     820             :   GEN a, p;
     821             :   long i, v;
     822             : 
     823      143024 :   switch(typ(x))
     824             :   {
     825             :     case t_INT:
     826      142716 :       return Z_issquare(x);
     827             : 
     828             :     case t_REAL:
     829          14 :       return (signe(x)>=0);
     830             : 
     831             :     case t_INTMOD:
     832          70 :       return Zn_ispower(gel(x,2), gel(x,1), gen_2, NULL);
     833             : 
     834             :     case t_FRAC:
     835          21 :       return Z_issquare(gel(x,1)) && Z_issquare(gel(x,2));
     836             : 
     837           7 :     case t_FFELT: return FF_issquareall(x, NULL);
     838             : 
     839             :     case t_COMPLEX:
     840           7 :       return 1;
     841             : 
     842             :     case t_PADIC:
     843          77 :       a = gel(x,4); if (!signe(a)) return 1;
     844          77 :       if (valp(x)&1) return 0;
     845          63 :       p = gel(x,2);
     846          63 :       if (!equaliu(p, 2)) return (kronecker(a,p) != -1);
     847             : 
     848          21 :       v = precp(x); /* here p=2, a is odd */
     849          21 :       if ((v>=3 && mod8(a) != 1 ) ||
     850           0 :           (v==2 && mod4(a) != 1)) return 0;
     851           7 :       return 1;
     852             : 
     853             :     case t_POLMOD:
     854          14 :       return polmodispower(x, gen_2, NULL);
     855             : 
     856             :     case t_POL:
     857          70 :       return polissquareall(x,NULL);
     858             : 
     859             :     case t_SER:
     860          21 :       if (!signe(x)) return 1;
     861          14 :       if (valp(x)&1) return 0;
     862           7 :       return issquare(gel(x,2));
     863             : 
     864             :     case t_RFRAC:
     865           7 :       av = avma; i = issquare(gmul(gel(x,1),gel(x,2)));
     866           7 :       avma = av; return i;
     867             :   }
     868           0 :   pari_err_TYPE("issquare",x);
     869           0 :   return 0; /* not reached */
     870             : }
     871           0 : GEN gissquare(GEN x) { return issquare(x)? gen_1: gen_0; }
     872           0 : GEN gissquareall(GEN x, GEN *pt) { return issquareall(x,pt)? gen_1: gen_0; }
     873             : 
     874             : long
     875        1386 : ispolygonal(GEN x, GEN S, GEN *N)
     876             : {
     877        1386 :   pari_sp av = avma;
     878             :   GEN D, d, n;
     879        1386 :   if (typ(x) != t_INT) pari_err_TYPE("ispolygonal", x);
     880        1386 :   if (typ(S) != t_INT) pari_err_TYPE("ispolygonal", S);
     881        1386 :   if (cmpiu(S,3) < 0) pari_err_DOMAIN("ispolygonal","s","<", utoipos(3),S);
     882        1386 :   if (signe(x) < 0) return 0;
     883        1386 :   if (signe(x) == 0) { if (N) *N = gen_0; return 1; }
     884        1260 :   if (is_pm1(x)) { if (N) *N = gen_1; return 1; }
     885             :   /* n = (sqrt( (8s - 16) x + (s-4)^2 ) + s - 4) / 2(s - 2) */
     886        1134 :   if (cmpiu(S, 1<<16) < 0) /* common case ! */
     887             :   {
     888         441 :     ulong s = S[2], r;
     889         504 :     if (s == 4) return Z_issquareall(x, N);
     890         378 :     if (s == 3)
     891           0 :       D = addiu(shifti(x, 3), 1);
     892             :     else
     893         378 :       D = addiu(mului(8*s - 16, x), (s-4)*(s-4));
     894         378 :     if (!Z_issquareall(D, &d)) { avma = av; return 0; }
     895         378 :     if (s == 3)
     896           0 :       d = subiu(d, 1);
     897             :     else
     898         378 :       d = addiu(d, s - 4);
     899         378 :     n = diviu_rem(d, 2*s - 4, &r);
     900         378 :     if (r) { avma = av; return 0; }
     901             :   }
     902             :   else
     903             :   {
     904         693 :     GEN r, S_2 = subiu(S,2), S_4 = subiu(S,4);
     905         693 :     D = addii(mulii(shifti(S_2,3), x), sqri(S_4));
     906         693 :     if (!Z_issquareall(D, &d)) { avma = av; return 0; }
     907         693 :     d = addii(d, S_4);
     908         693 :     n = dvmdii(shifti(d,-1), S_2, &r);
     909         693 :     if (r != gen_0) { avma = av; return 0; }
     910             :   }
     911        1071 :   if (N) *N = gerepileuptoint(av, n); else avma = av;
     912        1071 :   return 1;
     913             : }
     914             : 
     915             : /*********************************************************************/
     916             : /**                                                                 **/
     917             : /**                        PERFECT POWER                            **/
     918             : /**                                                                 **/
     919             : /*********************************************************************/
     920             : static long
     921         497 : polispower(GEN x, GEN K, GEN *pt)
     922             : {
     923             :   pari_sp av;
     924         497 :   long v, d, k = itos(K);
     925             :   GEN y, a, b;
     926             : 
     927         497 :   if (!signe(x))
     928             :   {
     929           7 :     if (pt) *pt = gcopy(x);
     930           7 :     return 1;
     931             :   }
     932         490 :   if (degpol(x) % k) return 0; /* degree not multiple of k */
     933         483 :   av = avma;
     934         483 :   y = NULL; /*-Wall*/
     935         483 :   v = RgX_valrem(x, &x);
     936         483 :   if (v % k) return 0;
     937         476 :   v /= k;
     938         476 :   a = gel(x,2); b = NULL;
     939         476 :   if (!ispower(a, K, &b)) { avma = av; return 0; }
     940         462 :   d = degpol(x);
     941         462 :   if (d)
     942             :   {
     943         322 :     GEN p = characteristic(x);
     944         322 :     a = leading_coeff(x);
     945         350 :     if (!ispower(a, K, &b)) { avma = av; return 0; }
     946         322 :     x = RgX_normalize(x);
     947         322 :     if (signe(p))
     948             :     {
     949         168 :       GEN T0, T = NULL;
     950         168 :       if (!BPSW_isprime(p))
     951           0 :         pari_err_IMPL("ispower in non-prime characteristic");
     952         168 :       if (RgX_is_FpXQX(x,&T,&p))
     953             :       { /* over Fq */
     954         168 :         T0 = T;
     955         168 :         if (T && typ(T) == t_FFELT) T = FF_mod(T);
     956         168 :         x = RgX_to_FqX(x,T,p);
     957         168 :         if (!FqX_ispower(x, k, T,p, pt)) { avma = av; return 0; }
     958         140 :         if (pt)
     959             :         {
     960         140 :           y = *pt;
     961         140 :           if (!T) y = FpX_to_mod(y, p);
     962         105 :           else if (typ(T0) == t_FFELT)
     963          70 :             y = FqX_to_FFX(y, T0);
     964             :           else
     965             :           {
     966          35 :             T = FpX_to_mod(T, p);
     967          35 :             y = gmul(y, gmodulsg(1,T));
     968             :           }
     969             :         }
     970         140 :         goto END;
     971             :       }
     972           0 :       if (cmpii(p,K) <= 0)
     973           0 :         pari_err_IMPL("ispower(general t_POL) in small characteristic");
     974             :     }
     975         154 :     y = gtrunc(gsqrtn(RgX_to_ser(x,lg(x)), K, NULL, 0));
     976         154 :     if (!RgX_equal(powgi(y, K), x)) { avma = av; return 0; }
     977             :   }
     978             :   else
     979         140 :     y = pol_1(varn(x));
     980             : END:
     981         434 :   if (pt)
     982             :   {
     983         434 :     if (!gequal1(a))
     984             :     {
     985          14 :       if (!b) b = gsqrtn(a, K, NULL, DEFAULTPREC);
     986          14 :       y = gmul(b,y);
     987             :     }
     988         434 :     if (v) y = RgX_shift_shallow(y, v);
     989         434 :     *pt = gerepilecopy(av, y);
     990             :   }
     991           0 :   else avma = av;
     992         434 :   return 1;
     993             : }
     994             : 
     995             : long
     996        1512 : Z_ispowerall(GEN x, ulong k, GEN *pt)
     997             : {
     998        1512 :   long s = signe(x);
     999             :   ulong mask;
    1000        1512 :   if (!s) { if (pt) *pt = gen_0; return 1; }
    1001        1512 :   if (s > 0) {
    1002        1358 :     if (k == 2) return Z_issquareall(x, pt);
    1003         987 :     if (k == 3) { mask = 1; return !!is_357_power(x, pt, &mask); }
    1004         217 :     if (k == 5) { mask = 2; return !!is_357_power(x, pt, &mask); }
    1005         203 :     if (k == 7) { mask = 4; return !!is_357_power(x, pt, &mask); }
    1006         196 :     return is_kth_power(x, k, pt);
    1007             :   }
    1008         154 :   if (!odd(k)) return 0;
    1009         140 :   if (Z_ispowerall(absi(x), k, pt))
    1010             :   {
    1011         140 :     if (pt) *pt = negi(*pt);
    1012         140 :     return 1;
    1013             :   };
    1014           0 :   return 0;
    1015             : }
    1016             : 
    1017             : /* is x a K-th power mod p ? Assume p prime. */
    1018             : int
    1019         210 : Fp_ispower(GEN x, GEN K, GEN p)
    1020             : {
    1021         210 :   pari_sp av = avma;
    1022             :   GEN p_1;
    1023             :   long r;
    1024         210 :   x = modii(x, p);
    1025         210 :   if (!signe(x) || equali1(x)) { avma = av; return 1; }
    1026             :   /* implies p > 2 */
    1027          91 :   p_1 = subiu(p,1);
    1028          91 :   K = gcdii(K, p_1);
    1029          91 :   if (equaliu(K, 2)) { r = kronecker(x,p); avma = av; return (r > 0); }
    1030          35 :   x = Fp_pow(x, diviiexact(p_1,K), p);
    1031          35 :   avma = av; return equali1(x);
    1032             : }
    1033             : 
    1034             : /* x unit defined modulo 2^e, e > 0, p prime */
    1035             : static int
    1036        2373 : U2_issquare(GEN x, long e)
    1037             : {
    1038        2373 :   long r = signe(x)>=0?mod8(x):8-mod8(x);
    1039        2373 :   if (e==1) return 1;
    1040        2373 :   if (e==2) return (r&3L) == 1;
    1041        2009 :   return r == 1;
    1042             : }
    1043             : /* x unit defined modulo p^e, e > 0, p prime */
    1044             : static int
    1045        4690 : Up_issquare(GEN x, GEN p, long e)
    1046        4690 : { return (equaliu(p,2))? U2_issquare(x, e): kronecker(x,p)==1; }
    1047             : 
    1048             : long
    1049        2548 : Zn_issquare(GEN d, GEN fn)
    1050             : {
    1051             :   long j, np;
    1052        2548 :   if (typ(d) != t_INT) pari_err_TYPE("Zn_issquare",d);
    1053        2548 :   if (typ(fn) == t_INT) return Zn_ispower(d, fn, gen_2, NULL);
    1054             :   /* integer factorization */
    1055        2548 :   np = nbrows(fn);
    1056        5320 :   for (j = 1; j <= np; ++j)
    1057             :   {
    1058        4970 :     GEN  r, p = gcoeff(fn, j, 1);
    1059        4970 :     long e = itos(gcoeff(fn, j, 2));
    1060        4970 :     long v = Z_pvalrem(d,p,&r);
    1061        4970 :     if (v < e && (odd(v) || !Up_issquare(r, p, e-v))) return 0;
    1062             :   }
    1063         350 :   return 1;
    1064             : }
    1065             : 
    1066             : static long
    1067        1113 : Qp_ispower(GEN x, GEN K, GEN *pt)
    1068             : {
    1069        1113 :   pari_sp av = avma;
    1070        1113 :   GEN z = Qp_sqrtn(x, K, NULL);
    1071        1113 :   if (!z) { avma = av; return 0; }
    1072         819 :   if (pt) *pt = z;
    1073         819 :   return 1;
    1074             : }
    1075             : 
    1076             : long
    1077     7002961 : ispower(GEN x, GEN K, GEN *pt)
    1078             : {
    1079             :   GEN z;
    1080             : 
    1081     7002961 :   if (!K) return gisanypower(x, pt);
    1082        2800 :   if (typ(K) != t_INT) pari_err_TYPE("ispower",K);
    1083        2800 :   if (signe(K) <= 0) pari_err_DOMAIN("ispower","exponent","<=",gen_0,K);
    1084        2800 :   if (equali1(K)) { if (pt) *pt = gcopy(x); return 1; }
    1085        2772 :   switch(typ(x)) {
    1086             :     case t_INT:
    1087         728 :       return Z_ispowerall(x, itou(K), pt);
    1088             :     case t_FRAC:
    1089             :     {
    1090          21 :       GEN a = gel(x,1), b = gel(x,2);
    1091          21 :       ulong k = itou(K);
    1092          21 :       if (pt) {
    1093          14 :         z = cgetg(3, t_FRAC);
    1094          14 :         if (Z_ispowerall(a, k, &a) && Z_ispowerall(b, k, &b)) {
    1095          14 :           *pt = z; gel(z,1) = a; gel(z,2) = b; return 1;
    1096             :         }
    1097           0 :         avma = (pari_sp)(z + 3); return 0;
    1098             :       }
    1099           7 :       return Z_ispower(a, k) && Z_ispower(b, k);
    1100             :     }
    1101             :     case t_INTMOD:
    1102         189 :       return Zn_ispower(gel(x,2), gel(x,1), K, pt);
    1103             :     case t_FFELT:
    1104         112 :       return FF_ispower(x, K, pt);
    1105             : 
    1106             :     case t_PADIC:
    1107        1113 :       return Qp_ispower(x, K, pt);
    1108             :     case t_POLMOD:
    1109          91 :       return polmodispower(x, K, pt);
    1110             :     case t_POL:
    1111         490 :       return polispower(x, K, pt);
    1112             :     case t_RFRAC: {
    1113           7 :       GEN a = gel(x,1), b = gel(x,2);
    1114           7 :       if (pt) {
    1115           7 :         z = cgetg(3, t_RFRAC);
    1116           7 :         if (ispower(a, K, &a) && polispower(b, K, &b)) {
    1117           7 :           *pt = z; gel(z,1) = a; gel(z,2) = b; return 1;
    1118             :         }
    1119           0 :         avma = (pari_sp)(z + 3); return 0;
    1120             :       }
    1121           0 :       return (ispower(a, K, NULL) && polispower(b, K, NULL));
    1122             :     }
    1123             :     case t_REAL:
    1124           7 :       if (signe(x) < 0 && !mpodd(K)) return 0;
    1125             :     case t_COMPLEX:
    1126          14 :       if (pt) *pt = gsqrtn(x, K, NULL, DEFAULTPREC);
    1127          14 :       return 1;
    1128             : 
    1129             :     case t_SER:
    1130           7 :       if (signe(x) && (!dvdsi(valp(x), K) || !ispower(gel(x,2), K, NULL)))
    1131           0 :         return 0;
    1132           7 :       if (pt) *pt = gsqrtn(x, K, NULL, DEFAULTPREC);
    1133           7 :       return 1;
    1134             :   }
    1135           0 :   pari_err_TYPE("ispower",x);
    1136           0 :   return 0; /* not reached */
    1137             : }
    1138             : 
    1139             : long
    1140     7000161 : gisanypower(GEN x, GEN *pty)
    1141             : {
    1142     7000161 :   long tx = typ(x);
    1143             :   ulong k, h;
    1144     7000161 :   if (tx == t_INT) return Z_isanypower(x, pty);
    1145          14 :   if (tx == t_FRAC)
    1146             :   {
    1147          14 :     pari_sp av = avma;
    1148          14 :     GEN fa, P, E, a = gel(x,1), b = gel(x,2);
    1149             :     long i, j, p, e;
    1150          14 :     int sw = (absi_cmp(a, b) > 0);
    1151             : 
    1152          14 :     if (sw) swap(a, b);
    1153          14 :     k = Z_isanypower(a, pty? &a: NULL);
    1154          14 :     if (!k)
    1155             :     { /* a = -1,1 or not a pure power */
    1156           7 :       if (!is_pm1(a)) { avma = av; return 0; }
    1157           7 :       if (signe(a) < 0) b = negi(b);
    1158           7 :       k = Z_isanypower(b, pty? &b: NULL);
    1159           7 :       if (!k || !pty) { avma = av; return k; }
    1160           7 :       *pty = gerepileupto(av, ginv(b));
    1161           7 :       return k;
    1162             :     }
    1163           7 :     fa = factoru(k);
    1164           7 :     P = gel(fa,1);
    1165           7 :     E = gel(fa,2); h = k;
    1166          14 :     for (i = lg(P) - 1; i > 0; i--)
    1167             :     {
    1168           7 :       p = P[i];
    1169           7 :       e = E[i];
    1170          21 :       for (j = 0; j < e; j++)
    1171          14 :         if (!is_kth_power(b, p, &b)) break;
    1172           7 :       if (j < e) k /= upowuu(p, e - j);
    1173             :     }
    1174           7 :     if (k == 1) { avma = av; return 0; }
    1175           7 :     if (!pty) { avma = av; return k; }
    1176           0 :     if (k != h) a = powiu(a, h/k);
    1177           0 :     *pty = gerepilecopy(av, mkfrac(a, b));
    1178           0 :     return k;
    1179             :   }
    1180           0 :   pari_err_TYPE("gisanypower", x);
    1181           0 :   return 0; /* not reached */
    1182             : }
    1183             : 
    1184             : /* v_p(x) = e != 0 for some p; return ispower(x,,&x), updating x.
    1185             :  * No need to optimize for 2,3,5,7 powers (done before) */
    1186             : static long
    1187      505729 : split_exponent(ulong e, GEN *x)
    1188             : {
    1189             :   GEN fa, P, E;
    1190      505729 :   long i, j, l, k = 1;
    1191      505729 :   if (e == 1) return 1;
    1192          14 :   fa = factoru(e);
    1193          14 :   P = gel(fa,1);
    1194          14 :   E = gel(fa,2); l = lg(P);
    1195          28 :   for (i = 1; i < l; i++)
    1196             :   {
    1197          14 :     ulong p = P[i];
    1198          28 :     for (j = 0; j < E[i]; j++)
    1199             :     {
    1200             :       GEN y;
    1201          14 :       if (!is_kth_power(*x, p, &y)) break;
    1202          14 :       k *= p; *x = y;
    1203             :     }
    1204             :   }
    1205          14 :   return k;
    1206             : }
    1207             : 
    1208             : static long
    1209      864316 : Z_isanypower_nosmalldiv(GEN *px)
    1210             : { /* any prime divisor of x is > 102 */
    1211      864316 :   const double LOG2_103 = 6.6865; /* lower bound for log_2(103) */
    1212      864316 :   const double LOG103 = 4.6347; /* lower bound for log(103) */
    1213             :   forprime_t T;
    1214      864316 :   ulong mask = 7, e2;
    1215             :   long k, ex;
    1216      864316 :   GEN y, x = *px;
    1217             : 
    1218      864316 :   k = 1;
    1219      864316 :   while (Z_issquareall(x, &y)) { k <<= 1; x = y; }
    1220      864316 :   while ( (ex = is_357_power(x, &y, &mask)) ) { k *= ex; x = y; }
    1221      864316 :   e2 = (ulong)((expi(x) + 1) / LOG2_103); /* >= log_103 (x) */
    1222      864316 :   if (u_forprime_init(&T, 11, e2))
    1223             :   {
    1224       16814 :     GEN logx = NULL;
    1225       16814 :     const ulong Q = 30011; /* prime */
    1226             :     ulong p, xmodQ;
    1227       16814 :     double dlogx = 0;
    1228             :     /* cut off at x^(1/p) ~ 2^30 bits which seems to be about optimum;
    1229             :      * for large p the modular checks are no longer competitively fast */
    1230       33649 :     while ( (ex = is_pth_power(x, &y, &T, 30)) )
    1231             :     {
    1232          21 :       k *= ex; x = y;
    1233          21 :       e2 = (ulong)((expi(x) + 1) / LOG2_103);
    1234          21 :       u_forprime_restrict(&T, e2);
    1235             :     }
    1236       16814 :     if (DEBUGLEVEL>4) err_printf("Z_isanypower: now k=%ld, x=%ld-bit\n", k, expi(x));
    1237       16814 :     xmodQ = umodiu(x, Q);
    1238             :     /* test Q | x, just in case */
    1239       16814 :     if (!xmodQ) { *px = x; return k * split_exponent(Z_lval(x,Q), px); }
    1240             :     /* x^(1/p) < 2^31 */
    1241       16800 :     p = T.p;
    1242       16800 :     if (p <= e2)
    1243             :     {
    1244       16786 :       logx = logr_abs( itor(x, DEFAULTPREC) );
    1245       16786 :       dlogx = rtodbl(logx);
    1246       16786 :       e2 = (ulong)(dlogx / LOG103); /* >= log_103(x) */
    1247             :     }
    1248      151480 :     while (p && p <= e2)
    1249             :     { /* is x a p-th power ? By computing y = round(x^(1/p)).
    1250             :        * Check whether y^p = x, first mod Q, then exactly. */
    1251      117880 :       pari_sp av = avma;
    1252             :       long e;
    1253      117880 :       GEN logy = divru(logx, p), y = grndtoi(mpexp(logy), &e);
    1254      117880 :       ulong ymodQ = umodiu(y,Q);
    1255      117880 :       if (e >= -10 || Fl_powu(ymodQ, p % (Q-1), Q) != xmodQ
    1256          21 :                    || !equalii(powiu(y, p), x)) avma = av;
    1257             :       else
    1258             :       {
    1259          21 :         k *= p; x = y; xmodQ = ymodQ; logx = logy; dlogx /= p;
    1260          21 :         e2 = (ulong)(dlogx / LOG103); /* >= log_103(x) */
    1261          21 :         u_forprime_restrict(&T, e2);
    1262          21 :         continue; /* if success, retry same p */
    1263             :       }
    1264      117859 :       p = u_forprime_next(&T);
    1265             :     }
    1266             :   }
    1267      864302 :   *px = x; return k;
    1268             : }
    1269             : 
    1270             : static ulong tinyprimes[] = {
    1271             :   2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71,
    1272             :   73, 79, 83, 89, 97, 101, 103, 107, 109, 113, 127, 131, 137, 139, 149, 151,
    1273             :   157, 163, 167, 173, 179, 181, 191, 193, 197, 199
    1274             : };
    1275             : 
    1276             : /* disregard the sign of x, caller will take care of x < 0 */
    1277             : static long
    1278     7000630 : Z_isanypower_aux(GEN x, GEN *pty)
    1279             : {
    1280             :   long ex, v, i, l, k;
    1281             :   GEN y, P, E;
    1282     7000630 :   ulong mask, e = 0;
    1283             : 
    1284     7000630 :   if (absi_cmp(x, gen_2) < 0) return 0; /* -1,0,1 */
    1285             : 
    1286     7000616 :   if (signe(x) < 0) x = negi(x);
    1287     7000616 :   k = l = 1;
    1288     7000616 :   P = cgetg(26 + 1, t_VECSMALL);
    1289     7000616 :   E = cgetg(26 + 1, t_VECSMALL);
    1290             :   /* trial division */
    1291   122926272 :   for(i = 0; i < 26; i++)
    1292             :   {
    1293    60128943 :     ulong p = tinyprimes[i];
    1294             :     int stop;
    1295    60128943 :     v = Z_lvalrem_stop(&x, p, &stop);
    1296    60128943 :     if (v)
    1297             :     {
    1298     7922439 :       P[l] = p;
    1299     7922439 :       E[l] = v; l++;
    1300    13588862 :       e = ugcd(e, v); if (e == 1) goto END;
    1301             :     }
    1302    54710628 :     if (stop) {
    1303      248108 :       if (is_pm1(x)) k = e;
    1304      248108 :       goto END;
    1305             :     }
    1306             :   }
    1307             : 
    1308     1334193 :   if (e)
    1309             :   { /* Bingo. Result divides e */
    1310             :     long v3, v5, v7;
    1311      505715 :     ulong e2 = e;
    1312      505715 :     v = u_lvalrem(e2, 2, &e2);
    1313      505715 :     if (v)
    1314             :     {
    1315      375277 :       for (i = 0; i < v; i++)
    1316             :       {
    1317      374185 :         if (!Z_issquareall(x, &y)) break;
    1318        1302 :         k <<= 1; x = y;
    1319             :       }
    1320             :     }
    1321      505715 :     mask = 0;
    1322      505715 :     v3 = u_lvalrem(e2, 3, &e2); if (v3) mask = 1;
    1323      505715 :     v5 = u_lvalrem(e2, 5, &e2); if (v5) mask |= 2;
    1324      505715 :     v7 = u_lvalrem(e2, 7, &e2); if (v7) mask |= 4;
    1325     1011507 :     while ( (ex = is_357_power(x, &y, &mask)) ) {
    1326          77 :       x = y;
    1327          77 :       switch(ex)
    1328             :       {
    1329          28 :         case 3: k *= 3; if (--v3 == 0) mask &= ~1; break;
    1330          28 :         case 5: k *= 5; if (--v5 == 0) mask &= ~2; break;
    1331          21 :         case 7: k *= 7; if (--v7 == 0) mask &= ~4; break;
    1332             :       }
    1333             :     }
    1334      505715 :     k *= split_exponent(e2, &x);
    1335             :   }
    1336             :   else
    1337      828478 :     k = Z_isanypower_nosmalldiv(&x);
    1338             : END:
    1339     7000616 :   if (pty && k != 1)
    1340             :   {
    1341        7959 :     if (e)
    1342             :     { /* add missing small factors */
    1343        6881 :       y = powuu(P[1], E[1] / k);
    1344        6881 :       for (i = 2; i < l; i++) y = mulii(y, powuu(P[i], E[i] / k));
    1345        6881 :       x = equali1(x)? y: mulii(x,y);
    1346             :     }
    1347        7959 :     *pty = x;
    1348             :   }
    1349     7000616 :   return k == 1? 0: k;
    1350             : }
    1351             : 
    1352             : long
    1353     7000630 : Z_isanypower(GEN x, GEN *pty)
    1354             : {
    1355     7000630 :   pari_sp av = avma;
    1356     7000630 :   long k = Z_isanypower_aux(x, pty);
    1357     7000630 :   if (!k) { avma = av; return 0; }
    1358        8001 :   if (signe(x) < 0)
    1359             :   {
    1360          42 :     long v = vals(k);
    1361          42 :     if (v)
    1362             :     {
    1363             :       GEN y;
    1364          28 :       k >>= v;
    1365          28 :       if (k == 1) { avma = av; return 0; }
    1366          21 :       if (!pty) { avma = av; return k; }
    1367          14 :       y = *pty;
    1368          14 :       y = powiu(y, 1<<v);
    1369          14 :       togglesign(y);
    1370          14 :       *pty = gerepileuptoint(av, y);
    1371          14 :       return k;
    1372             :     }
    1373          14 :     if (pty) togglesign_safe(pty);
    1374             :   }
    1375        7973 :   if (pty) *pty = gerepilecopy(av, *pty); else avma = av;
    1376        7973 :   return k;
    1377             : }
    1378             : 
    1379             : /* Faster than */
    1380             : /*   !cmpii(n, int2n(vali(n))) */
    1381             : /*   !cmpis(shifti(n, -vali(n)), 1) */
    1382             : /*   expi(n) == vali(n) */
    1383             : /*   hamming(n) == 1 */
    1384             : /* even for single-word values, and much faster for multiword values. */
    1385             : /* If all you have is a word, you can just use n & !(n & (n-1)). */
    1386             : long
    1387       30824 : Z_ispow2(GEN n)
    1388             : {
    1389             :   GEN xp;
    1390             :   long i, lx;
    1391             :   ulong u;
    1392       30824 :   if (signe(n) != 1) return 0;
    1393       30817 :   xp = int_LSW(n);
    1394       30817 :   lx = lgefint(n);
    1395       30817 :   u = *xp;
    1396       30849 :   for (i = 3; i < lx; ++i)
    1397             :   {
    1398       28366 :     if (u) return 0;
    1399          32 :     xp = int_nextW(xp);
    1400          32 :     u = *xp;
    1401             :   }
    1402        2483 :   return !(u & (u-1)); /* faster than hamming_word(u) == 1 */
    1403             : }
    1404             : 
    1405             : static long
    1406      841627 : isprimepower_i(GEN n, GEN *pt, long flag)
    1407             : {
    1408      841627 :   pari_sp av = avma;
    1409             :   long i, v;
    1410             : 
    1411      841627 :   if (typ(n) != t_INT) pari_err_TYPE("isprimepower", n);
    1412      841626 :   if (signe(n) <= 0) return 0;
    1413             : 
    1414      841626 :   if (lgefint(n) == 3)
    1415             :   {
    1416             :     ulong p;
    1417      541061 :     v = uisprimepower(n[2], &p);
    1418      541061 :     if (v)
    1419             :     {
    1420       54849 :       if (pt) *pt = utoipos(p);
    1421       54849 :       return v;
    1422             :     }
    1423      486212 :     return 0;
    1424             :   }
    1425     1662025 :   for (i = 0; i < 26; i++)
    1426             :   {
    1427     1626187 :     ulong p = tinyprimes[i];
    1428     1626187 :     v = Z_lvalrem(n, p, &n);
    1429     1626188 :     if (v)
    1430             :     {
    1431      264728 :       avma = av;
    1432      264728 :       if (!is_pm1(n)) return 0;
    1433         332 :       if (pt) *pt = utoipos(p);
    1434         331 :       return v;
    1435             :     }
    1436             :   }
    1437             :   /* p | n => p >= 103 */
    1438       35838 :   v = Z_isanypower_nosmalldiv(&n); /* expensive */
    1439       35838 :   if (!(flag? isprime(n): BPSW_psp(n))) { avma = av; return 0; }
    1440        5526 :   if (pt) *pt = gerepilecopy(av, n); else avma = av;
    1441        5526 :   return v;
    1442             : }
    1443             : long
    1444      840098 : isprimepower(GEN n, GEN *pt) { return isprimepower_i(n,pt,1); }
    1445             : long
    1446        1529 : ispseudoprimepower(GEN n, GEN *pt) { return isprimepower_i(n,pt,0); }
    1447             : 
    1448             : long
    1449      541670 : uisprimepower(ulong n, ulong *pp)
    1450             : { /* We must have CUTOFF^11 >= ULONG_MAX and CUTOFF^3 < ULONG_MAX.
    1451             :    * Tests suggest that 200-300 is the best range for 64-bit platforms. */
    1452      541670 :   const ulong CUTOFF = 200UL;
    1453      541670 :   const long TINYCUTOFF = 46;  /* tinyprimes[45] = 199 */
    1454      541670 :   const ulong CUTOFF3 = CUTOFF*CUTOFF*CUTOFF;
    1455             : #ifdef LONG_IS_64BIT
    1456             :   /* primes preceeding the appropriate root of ULONG_MAX. */
    1457      481440 :   const ulong ROOT9 = 137;
    1458      481440 :   const ulong ROOT8 = 251;
    1459      481440 :   const ulong ROOT7 = 563;
    1460      481440 :   const ulong ROOT5 = 7129;
    1461      481440 :   const ulong ROOT4 = 65521;
    1462             : #else
    1463       60230 :   const ulong ROOT9 = 11;
    1464       60230 :   const ulong ROOT8 = 13;
    1465       60230 :   const ulong ROOT7 = 23;
    1466       60230 :   const ulong ROOT5 = 83;
    1467       60230 :   const ulong ROOT4 = 251;
    1468             : #endif
    1469             :   ulong mask;
    1470             :   long v, i;
    1471             :   int e;
    1472      541670 :   if (n < 2) return 0;
    1473      541656 :   if (!odd(n)) {
    1474      270683 :     if (n & (n-1)) return 0;
    1475         886 :     *pp = 2; return vals(n);
    1476             :   }
    1477      270973 :   if (n < 8) { *pp = n; return 1; } /* 3,5,7 */
    1478     3653596 :   for (i = 1/*skip p=2*/; i < TINYCUTOFF; i++)
    1479             :   {
    1480     3594527 :     ulong p = tinyprimes[i];
    1481     3594527 :     if (n % p == 0)
    1482             :     {
    1483      211519 :       v = u_lvalrem(n, p, &n);
    1484      211519 :       if (n == 1) { *pp = p; return v; }
    1485      209438 :       return 0;
    1486             :     }
    1487             :   }
    1488             :   /* p | n => p >= CUTOFF */
    1489             : 
    1490       59069 :   if (n < CUTOFF3)
    1491             :   {
    1492       46354 :     if (n < CUTOFF*CUTOFF || uisprime_101(n)) { *pp = n; return 1; }
    1493           0 :     if (uissquareall(n, &n)) { *pp = n; return 2; }
    1494           0 :     return 0;
    1495             :   }
    1496             : 
    1497             :   /* Check for squares, fourth powers, and eighth powers as appropriate. */
    1498       12715 :   v = 1;
    1499       12715 :   if (uissquareall(n, &n)) {
    1500           0 :     v <<= 1;
    1501           0 :     if (CUTOFF <= ROOT4 && uissquareall(n, &n)) {
    1502           0 :       v <<= 1;
    1503           0 :       if (CUTOFF <= ROOT8 && uissquareall(n, &n)) v <<= 1;
    1504             :     }
    1505             :   }
    1506             : 
    1507       12715 :   if (CUTOFF > ROOT5) mask = 1;
    1508             :   else
    1509             :   {
    1510       12714 :     const ulong CUTOFF5 = CUTOFF3*CUTOFF*CUTOFF;
    1511       12714 :     if (n < CUTOFF5) mask = 1; else mask = 3;
    1512       12714 :     if (CUTOFF <= ROOT7)
    1513             :     {
    1514       12714 :       const ulong CUTOFF7 = CUTOFF5*CUTOFF*CUTOFF;
    1515       12714 :       if (n >= CUTOFF7) mask = 7;
    1516             :     }
    1517             :   }
    1518             : 
    1519       12715 :   if (CUTOFF <= ROOT9 && (e = uis_357_power(n, &n, &mask))) { v *= e; mask=1; }
    1520       12715 :   if ((e = uis_357_power(n, &n, &mask))) v *= e;
    1521             : 
    1522       12715 :   if (uisprime_101(n)) { *pp = n; return v; }
    1523        6984 :   return 0;
    1524             : }
    1525             : 
    1526             : /*********************************************************************/
    1527             : /**                                                                 **/
    1528             : /**                        KRONECKER SYMBOL                         **/
    1529             : /**                                                                 **/
    1530             : /*********************************************************************/
    1531             : /* t = 3,5 mod 8 ?  (= 2 not a square mod t) */
    1532             : static int
    1533   539671898 : ome(long t)
    1534             : {
    1535   539671898 :   switch(t & 7)
    1536             :   {
    1537             :     case 3:
    1538   310357593 :     case 5: return 1;
    1539   229314305 :     default: return 0;
    1540             :   }
    1541             : }
    1542             : /* t a t_INT, is t = 3,5 mod 8 ? */
    1543             : static int
    1544     4412282 : gome(GEN t)
    1545     4412282 : { return signe(t)? ome( mod2BIL(t) ): 0; }
    1546             : 
    1547             : /* t a t_INT, return 1 if t = 3 (mod 4), 0 otherwise */
    1548             : static int
    1549       18445 : eps(GEN t)
    1550             : {
    1551       18445 :   switch(signe(t))
    1552             :   {
    1553        4977 :     case -1: return mod4(t) == 1;
    1554       13468 :     case 1:  return mod4(t) == 3;
    1555           0 :     default: return 0;
    1556             :   }
    1557             : }
    1558             : 
    1559             : /* assume y odd, return kronecker(x,y) * s */
    1560             : static long
    1561   407723261 : krouu_s(ulong x, ulong y, long s)
    1562             : {
    1563   407723261 :   ulong x1 = x, y1 = y, z;
    1564  2196968027 :   while (x1)
    1565             :   {
    1566  1381610123 :     long r = vals(x1);
    1567  1381704189 :     if (r)
    1568             :     {
    1569   755234115 :       if (odd(r) && ome(y1)) s = -s;
    1570   755051431 :       x1 >>= r;
    1571             :     }
    1572  1381521505 :     if (x1 & y1 & 2) s = -s;
    1573  1381521505 :     z = y1 % x1; y1 = x1; x1 = z;
    1574             :   }
    1575   407634643 :   return (y1 == 1)? s: 0;
    1576             : }
    1577             : 
    1578             : long
    1579     5054204 : kronecker(GEN x, GEN y)
    1580             : {
    1581     5054204 :   pari_sp av = avma;
    1582     5054204 :   long s = 1, r;
    1583             :   ulong xu, yu;
    1584             : 
    1585     5054204 :   if (typ(x) != t_INT) pari_err_TYPE("kronecker",x);
    1586     5054204 :   if (typ(y) != t_INT) pari_err_TYPE("kronecker",y);
    1587     5054204 :   switch (signe(y))
    1588             :   {
    1589           0 :     case -1: y = negi(y); if (signe(x) < 0) s = -1; break;
    1590           0 :     case 0: return is_pm1(x);
    1591             :   }
    1592     5054204 :   r = vali(y);
    1593     5054204 :   if (r)
    1594             :   {
    1595       10403 :     if (!mpodd(x)) { avma = av; return 0; }
    1596       10200 :     if (odd(r) && gome(x)) s = -s;
    1597       10200 :     y = shifti(y,-r);
    1598             :   }
    1599     5054001 :   x = modii(x,y);
    1600    10183442 :   while (lgefint(x) > 3) /* x < y */
    1601             :   {
    1602             :     GEN z;
    1603       75440 :     r = vali(x);
    1604       75440 :     if (r)
    1605             :     {
    1606       41328 :       if (odd(r) && gome(y)) s = -s;
    1607       41328 :       x = shifti(x,-r);
    1608             :     }
    1609             :     /* x=3 mod 4 && y=3 mod 4 ? (both are odd here) */
    1610       75440 :     if (mod2BIL(x) & mod2BIL(y) & 2) s = -s;
    1611       75440 :     z = remii(y,x); y = x; x = z;
    1612       75440 :     if (gc_needed(av,2))
    1613             :     {
    1614           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"kronecker");
    1615           0 :       gerepileall(av, 2, &x, &y);
    1616             :     }
    1617             :   }
    1618     5054001 :   xu = itou(x);
    1619     5054001 :   if (!xu) return is_pm1(y)? s: 0;
    1620     5033336 :   r = vals(xu);
    1621     5033336 :   if (r)
    1622             :   {
    1623     3492835 :     if (odd(r) && gome(y)) s = -s;
    1624     3492835 :     xu >>= r;
    1625             :   }
    1626             :   /* x=3 mod 4 && y=3 mod 4 ? (both are odd here) */
    1627     5033336 :   if (xu & mod2BIL(y) & 2) s = -s;
    1628     5033336 :   yu = umodiu(y, xu);
    1629     5033336 :   avma = av; return krouu_s(yu, xu, s);
    1630             : }
    1631             : 
    1632             : long
    1633       27671 : krois(GEN x, long y)
    1634             : {
    1635             :   ulong yu;
    1636       27671 :   long s = 1;
    1637             : 
    1638       27671 :   if (y <= 0)
    1639             :   {
    1640           0 :     if (y == 0) return is_pm1(x);
    1641           0 :     yu = (ulong)-y; if (signe(x) < 0) s = -1;
    1642             :   }
    1643             :   else
    1644       27671 :     yu = (ulong)y;
    1645       27671 :   if (!odd(yu))
    1646             :   {
    1647             :     long r;
    1648       13335 :     if (!mpodd(x)) return 0;
    1649        9961 :     r = vals(yu); yu >>= r;
    1650        9961 :     if (odd(r) && gome(x)) s = -s;
    1651             :   }
    1652       24297 :   return krouu_s(umodiu(x, yu), yu, s);
    1653             : }
    1654             : /* assume y != 0 */
    1655             : long
    1656   283954419 : kroiu(GEN x, ulong y)
    1657             : {
    1658             :   long r;
    1659   283954419 :   if (odd(y)) return krouu_s(umodiu(x,y), y, 1);
    1660     1800022 :   if (!mpodd(x)) return 0;
    1661     1764042 :   r = vals(y); y >>= r;
    1662     1764042 :   return krouu_s(umodiu(x,y), y, (odd(r) && gome(x))? -1: 1);
    1663             : }
    1664             : 
    1665             : /* assume y > 0, odd, return s * kronecker(x,y) */
    1666             : static long
    1667       31274 : krouodd(ulong x, GEN y, long s)
    1668             : {
    1669             :   long r;
    1670       31274 :   if (lgefint(y) == 3) return krouu_s(x, y[2], s);
    1671       13046 :   if (!x) return 0; /* y != 1 */
    1672       13046 :   r = vals(x);
    1673       13046 :   if (r)
    1674             :   {
    1675        8429 :     if (odd(r) && gome(y)) s = -s;
    1676        8429 :     x >>= r;
    1677             :   }
    1678             :   /* x=3 mod 4 && y=3 mod 4 ? (both are odd here) */
    1679       13046 :   if (x & mod2BIL(y) & 2) s = -s;
    1680       13046 :   return krouu_s(umodiu(y,x), x, s);
    1681             : }
    1682             : 
    1683             : long
    1684       31162 : krosi(long x, GEN y)
    1685             : {
    1686       31162 :   const pari_sp av = avma;
    1687       31162 :   long s = 1, r;
    1688       31162 :   switch (signe(y))
    1689             :   {
    1690           0 :     case -1: y = negi(y); if (x < 0) s = -1; break;
    1691           0 :     case 0: return (x==1 || x==-1);
    1692             :   }
    1693       31162 :   r = vali(y);
    1694       31162 :   if (r)
    1695             :   {
    1696           0 :     if (!odd(x)) { avma = av; return 0; }
    1697           0 :     if (odd(r) && ome(x)) s = -s;
    1698           0 :     y = shifti(y,-r);
    1699             :   }
    1700       31162 :   if (x < 0) { x = -x; if (mod4(y) == 3) s = -s; }
    1701       31162 :   s = krouodd((ulong)x, y, s);
    1702       31162 :   avma = av; return s;
    1703             : }
    1704             : 
    1705             : long
    1706         112 : kroui(ulong x, GEN y)
    1707             : {
    1708         112 :   const pari_sp av = avma;
    1709         112 :   long s = 1, r;
    1710         112 :   switch (signe(y))
    1711             :   {
    1712           0 :     case -1: y = negi(y); break;
    1713           0 :     case 0: return x==1UL;
    1714             :   }
    1715         112 :   r = vali(y);
    1716         112 :   if (r)
    1717             :   {
    1718           0 :     if (!odd(x)) { avma = av; return 0; }
    1719           0 :     if (odd(r) && ome(x)) s = -s;
    1720           0 :     y = shifti(y,-r);
    1721             :   }
    1722         112 :   s = krouodd(x, y, s);
    1723         112 :   avma = av; return s;
    1724             : }
    1725             : 
    1726             : long
    1727    46125977 : kross(long x, long y)
    1728             : {
    1729             :   ulong yu;
    1730    46125977 :   long s = 1;
    1731             : 
    1732    46125977 :   if (y <= 0)
    1733             :   {
    1734         385 :     if (y == 0) return (labs(x)==1);
    1735         385 :     yu = (ulong)-y; if (x < 0) s = -1;
    1736             :   }
    1737             :   else
    1738    46125592 :     yu = (ulong)y;
    1739    46125977 :   if (!odd(yu))
    1740             :   {
    1741             :     long r;
    1742     4745297 :     if (!odd(x)) return 0;
    1743     4744618 :     r = vals(yu); yu >>= r;
    1744     4744618 :     if (odd(r) && ome(x)) s = -s;
    1745             :   }
    1746    46125298 :   x %= (long)yu; if (x < 0) x += yu;
    1747    46125298 :   return krouu_s((ulong)x, yu, s);
    1748             : }
    1749             : 
    1750             : long
    1751    72479673 : krouu(ulong x, ulong y)
    1752             : {
    1753             :   long r;
    1754    72479673 :   if (odd(y)) return krouu_s(x, y, 1);
    1755        1884 :   if (!odd(x)) return 0;
    1756        1884 :   r = vals(y); y >>= r;
    1757        1884 :   return krouu_s(x, y, (odd(r) && ome(x))? -1: 1);
    1758             : }
    1759             : 
    1760             : /*********************************************************************/
    1761             : /**                                                                 **/
    1762             : /**                          HILBERT SYMBOL                         **/
    1763             : /**                                                                 **/
    1764             : /*********************************************************************/
    1765             : /* x,y are t_INT or t_REAL */
    1766             : static long
    1767        9982 : mphilbertoo(GEN x, GEN y)
    1768             : {
    1769        9982 :   long sx = signe(x), sy = signe(y);
    1770        9982 :   if (!sx || !sy) return 0;
    1771        9982 :   return (sx < 0 && sy < 0)? -1: 1;
    1772             : }
    1773             : 
    1774             : long
    1775       53144 : hilbertii(GEN x, GEN y, GEN p)
    1776             : {
    1777             :   pari_sp av;
    1778             :   long oddvx, oddvy, z;
    1779             : 
    1780       53144 :   if (!p) return mphilbertoo(x,y);
    1781       43183 :   if (is_pm1(p) || signe(p) < 0) pari_err_PRIME("hilbertii",p);
    1782       43183 :   if (!signe(x) || !signe(y)) return 0;
    1783       43162 :   av = avma;
    1784       43162 :   oddvx = odd(Z_pvalrem(x,p,&x));
    1785       43162 :   oddvy = odd(Z_pvalrem(y,p,&y));
    1786             :   /* x, y are p-units, compute hilbert(x * p^oddvx, y * p^oddvy, p) */
    1787       43162 :   if (equaliu(p, 2))
    1788             :   {
    1789       10689 :     z = (eps(x) && eps(y))? -1: 1;
    1790       10689 :     if (oddvx && gome(y)) z = -z;
    1791       10689 :     if (oddvy && gome(x)) z = -z;
    1792             :   }
    1793             :   else
    1794             :   {
    1795       32473 :     z = (oddvx && oddvy && eps(p))? -1: 1;
    1796       32473 :     if (oddvx && kronecker(y,p) < 0) z = -z;
    1797       32473 :     if (oddvy && kronecker(x,p) < 0) z = -z;
    1798             :   }
    1799       43162 :   avma = av; return z;
    1800             : }
    1801             : 
    1802             : static void
    1803         196 : err_prec(void) { pari_err_PREC("hilbert"); }
    1804             : static void
    1805         161 : err_p(GEN p, GEN q) { pari_err_MODULUS("hilbert", p,q); }
    1806             : static void
    1807          56 : err_oo(GEN p) { pari_err_MODULUS("hilbert", p, strtoGENstr("oo")); }
    1808             : 
    1809             : /* x t_INTMOD, *pp = prime or NULL [ unset, set it to x.mod ].
    1810             :  * Return lift(x) provided it's p-adic accuracy is large enough to decide
    1811             :  * hilbert()'s value [ problem at p = 2 ] */
    1812             : static GEN
    1813         420 : lift_intmod(GEN x, GEN *pp)
    1814             : {
    1815         420 :   GEN p = *pp, N = gel(x,1);
    1816         420 :   x = gel(x,2);
    1817         420 :   if (!p)
    1818             :   {
    1819         266 :     *pp = p = N;
    1820         266 :     switch(itos_or_0(p))
    1821             :     {
    1822             :       case 2:
    1823         126 :       case 4: err_prec();
    1824             :     }
    1825         140 :     return x;
    1826             :   }
    1827         154 :   if (!signe(p)) err_oo(N);
    1828         112 :   if (equaliu(p,2))
    1829          42 :   { if (vali(N) <= 2) err_prec(); }
    1830             :   else
    1831          70 :   { if (!dvdii(N,p)) err_p(N,p); }
    1832          28 :   if (!signe(x)) err_prec();
    1833          21 :   return x;
    1834             : }
    1835             : /* x t_PADIC, *pp = prime or NULL [ unset, set it to x.p ].
    1836             :  * Return lift(x)*p^(v(x) mod 2) provided it's p-adic accuracy is large enough
    1837             :  * to decide hilbert()'s value [ problem at p = 2 ]*/
    1838             : static GEN
    1839         210 : lift_padic(GEN x, GEN *pp)
    1840             : {
    1841         210 :   GEN p = *pp, q = gel(x,2), y = gel(x,4);
    1842         210 :   if (!p) *pp = p = q;
    1843         147 :   else if (!equalii(p,q)) err_p(p, q);
    1844         105 :   if (equaliu(p,2) && precp(x) <= 2) err_prec();
    1845          70 :   if (!signe(y)) err_prec();
    1846          70 :   return odd(valp(x))? mulii(p,y): y;
    1847             : }
    1848             : 
    1849             : long
    1850         658 : hilbert(GEN x, GEN y, GEN p)
    1851             : {
    1852         658 :   pari_sp av = avma;
    1853         658 :   long tx = typ(x), ty = typ(y), z;
    1854             : 
    1855         658 :   if (p && typ(p) != t_INT) pari_err_TYPE("hilbert",p);
    1856         658 :   if (tx == t_REAL)
    1857             :   {
    1858          77 :     if (p && signe(p)) err_oo(p);
    1859          63 :     switch (ty)
    1860             :     {
    1861             :       case t_INT:
    1862           7 :       case t_REAL: return mphilbertoo(x,y);
    1863           0 :       case t_FRAC: return mphilbertoo(x,gel(y,1));
    1864          56 :       default: pari_err_TYPE2("hilbert",x,y);
    1865             :     }
    1866             :   }
    1867         581 :   if (ty == t_REAL)
    1868             :   {
    1869          14 :     if (p && signe(p)) err_oo(p);
    1870          14 :     switch (tx)
    1871             :     {
    1872             :       case t_INT:
    1873          14 :       case t_REAL: return mphilbertoo(x,y);
    1874           0 :       case t_FRAC: return mphilbertoo(gel(x,1),y);
    1875           0 :       default: pari_err_TYPE2("hilbert",x,y);
    1876             :     }
    1877             :   }
    1878         567 :   if (tx == t_INTMOD) { x = lift_intmod(x, &p); tx = t_INT; }
    1879         364 :   if (ty == t_INTMOD) { y = lift_intmod(y, &p); ty = t_INT; }
    1880             : 
    1881         308 :   if (tx == t_PADIC) { x = lift_padic(x, &p); tx = t_INT; }
    1882         245 :   if (ty == t_PADIC) { y = lift_padic(y, &p); ty = t_INT; }
    1883             : 
    1884         168 :   if (tx == t_FRAC) { tx = t_INT; x = p? mulii(gel(x,1),gel(x,2)): gel(x,1); }
    1885         168 :   if (ty == t_FRAC) { ty = t_INT; y = p? mulii(gel(y,1),gel(y,2)): gel(y,1); }
    1886             : 
    1887         168 :   if (tx != t_INT || ty != t_INT) pari_err_TYPE2("hilbert",x,y);
    1888         168 :   if (p && !signe(p)) p = NULL;
    1889         168 :   z = hilbertii(x,y,p); avma = av; return z;
    1890             : }
    1891             : 
    1892             : /*******************************************************************/
    1893             : /*                                                                 */
    1894             : /*                       SQUARE ROOT MODULO p                      */
    1895             : /*                                                                 */
    1896             : /*******************************************************************/
    1897             : 
    1898             : static ulong
    1899    22196598 : Fl_2gener_pre_all(long e, ulong p, ulong pi, ulong *pt_m)
    1900             : {
    1901             :   ulong y, m;
    1902             :   long k, i;
    1903    22196598 :   ulong q = (p-1) >> e; /* q = (p-1)/2^oo is odd */
    1904    39566725 :   for (k=2; ; k++)
    1905             :   { /* loop terminates for k < p (even if p composite) */
    1906    39566725 :     i = krouu(k, p);
    1907    39566725 :     if (i >= 0)
    1908             :     {
    1909    17370134 :       if (i) continue;
    1910           7 :       pari_err_PRIME("Fl_sqrt [modulus]",utoi(p));
    1911             :     }
    1912    22196591 :     y = m = Fl_powu_pre(k, q, p, pi);
    1913    60642158 :     for (i=1; i<e; i++)
    1914    38445567 :       if ((m = Fl_sqr_pre(m, p, pi)) == 1) break;
    1915    22196591 :     if (i == e) break; /* success */
    1916    17370127 :   }
    1917    22196591 :   *pt_m = m;
    1918    22196591 :   return y;
    1919             : }
    1920             : 
    1921             : /* Tonelli-Shanks. Assume p is prime and (a,p) != -1. */
    1922             : static ulong
    1923    53873656 : Fl_sqrt_i(ulong a, ulong p, ulong pi, ulong y, ulong m)
    1924             : {
    1925             :   long i, e, k;
    1926             :   ulong p1, q, v, w;
    1927             : 
    1928    53873656 :   if (!a) return 0;
    1929    53027525 :   p1 = p - 1; e = vals(p1);
    1930    53035824 :   if (e == 0) /* p = 2 */
    1931             :   {
    1932      417193 :     if (p != 2) pari_err_PRIME("Fl_sqrt [modulus]",utoi(p));
    1933      417186 :     return ((a & 1) == 0)? 0: 1;
    1934             :   }
    1935    52618631 :   q = p1 >> e; /* q = (p-1)/2^oo is odd */
    1936    52618631 :   if (e == 1)    y = p1;
    1937    22196598 :   else if (y==0) y = Fl_2gener_pre_all(e, p, pi, &m);
    1938    52618624 :   p1 = Fl_powu_pre(a, q >> 1, p, pi); /* a ^ [(q-1)/2] */
    1939    52607554 :   if (!p1) return 0;
    1940    52607554 :   v = Fl_mul_pre(a, p1, p, pi);
    1941    52612505 :   w = Fl_mul_pre(v, p1, p, pi);
    1942   123806399 :   while (w != 1)
    1943             :   { /* a*w = v^2, y primitive 2^e-th root of 1
    1944             :        a square --> w even power of y, hence w^(2^(e-1)) = 1 */
    1945    18679051 :     p1 = Fl_sqr_pre(w,p,pi);
    1946    18679053 :     for (k=1; p1 != 1 && k < e; k++) p1 = Fl_sqr_pre(p1,p,pi);
    1947    18679053 :     if (k == e) return ~0UL;
    1948             :     /* w ^ (2^k) = 1 --> w = y ^ (u * 2^(e-k)), u odd */
    1949    18583723 :     p1 = y;
    1950    18583723 :     for (i=1; i < e-k; i++) p1 = Fl_sqr_pre(p1, p, pi);
    1951    18583723 :     y = Fl_sqr_pre(p1, p, pi); e = k;
    1952    18583723 :     w = Fl_mul_pre(y, w, p, pi);
    1953    18583723 :     v = Fl_mul_pre(v, p1, p, pi);
    1954             :   }
    1955    52515597 :   p1 = p - v; if (v > p1) v = p1;
    1956    52515597 :   return v;
    1957             : }
    1958             : 
    1959             : ulong
    1960    50409050 : Fl_sqrt(ulong a, ulong p)
    1961             : {
    1962    50409050 :   ulong pi = get_Fl_red(p);
    1963    50402406 :   return Fl_sqrt_i(a, p, pi, 0, 0);
    1964             : }
    1965             : 
    1966             : ulong
    1967     3472431 : Fl_sqrt_pre(ulong a, ulong p, ulong pi)
    1968             : {
    1969     3472431 :   return Fl_sqrt_i(a, p, pi, 0, 0);
    1970             : }
    1971             : 
    1972             : static ulong
    1973       63230 : Fl_lgener_pre_all(ulong l, long e, ulong r, ulong p, ulong pi, ulong *pt_m)
    1974             : {
    1975             :   ulong x, y, m;
    1976       63230 :   ulong le1 = upowuu(l, e-1);
    1977      100463 :   for (x = 2; ; x++)
    1978             :   {
    1979      100463 :     y = Fl_powu_pre(x, r, p, pi);
    1980      100463 :     if (y==1) continue;
    1981       78196 :     m = Fl_powu_pre(y, le1, p, pi);
    1982       78196 :     if (m != 1) break;
    1983       37233 :   }
    1984       63230 :   *pt_m = m;
    1985       63230 :   return y;
    1986             : }
    1987             : 
    1988             : /* solve x^l = a , l prime in G of order q.
    1989             :  *
    1990             :  * q =  (l^e)*r, e >= 1, (r,l) = 1
    1991             :  * y generates the l-Sylow of G
    1992             :  * m = y^(l^(e-1)) != 1 */
    1993             : static ulong
    1994      115508 : Fl_sqrtl_i(ulong a, ulong l, ulong p, ulong pi, ulong y, ulong m)
    1995             : {
    1996             :   ulong p1, v, w, z, dl, zm;
    1997             :   ulong r, e, u2;
    1998      115508 :   if (a==0) return a;
    1999      115502 :   e = u_lvalrem(p-1, l, &r);
    2000      115501 :   u2 = Fl_inv(l%r, r);
    2001      115503 :   v = Fl_powu_pre(a, u2, p,pi);
    2002      115503 :   w = Fl_powu_pre(v, l, p,pi);
    2003      115503 :   w = Fl_mul_pre(w, Fl_inv(a, p),p,pi);
    2004      115503 :   if (w==1) return v;
    2005       63230 :   if (y==0) y = Fl_lgener_pre_all(l, e, r, p, pi, &m);
    2006      145426 :   while (w!=1)
    2007             :   {
    2008       67654 :     ulong k = 0;
    2009       67654 :     p1 = w;
    2010             :     do
    2011             :     {
    2012      101172 :       z = p1; p1 = Fl_powu_pre(p1, l, p, pi);
    2013      101172 :       k++;
    2014      101172 :     } while (p1!=1);
    2015       67654 :     if (k==e) return ~0UL;
    2016       18966 :     dl = 0; zm = 1;
    2017       66822 :     while (z!=zm)
    2018             :     {
    2019       28890 :       zm = Fl_mul_pre(zm, m, p, pi); dl++;
    2020             :     }
    2021       18966 :     dl = Fl_neg(dl, l);
    2022       18966 :     p1 = Fl_powu_pre(y,dl*upowuu(l,e-k-1),p,pi);
    2023       18966 :     m = Fl_powu_pre(m, dl, p, pi);
    2024       18966 :     e = k;
    2025       18966 :     v = Fl_mul_pre(p1,v,p,pi);
    2026       18966 :     y = Fl_powu_pre(p1,l,p,pi);
    2027       18966 :     w = Fl_mul_pre(y,w,p,pi);
    2028             :   }
    2029       14542 :   return v;
    2030             : }
    2031             : 
    2032             : ulong
    2033      115508 : Fl_sqrtl_pre(ulong a, ulong l, ulong p, ulong pi)
    2034             : {
    2035      115508 :   return Fl_sqrtl_i(a, l, p, pi, 0, 0);
    2036             : }
    2037             : 
    2038             : ulong
    2039           0 : Fl_sqrtl(ulong a, ulong l, ulong p)
    2040             : {
    2041           0 :   ulong pi = get_Fl_red(p);
    2042           0 :   return Fl_sqrtl_i(a, l, p, pi, 0, 0);
    2043             : }
    2044             : 
    2045             : /* Cipolla is better than Tonelli-Shanks when e = v_2(p-1) is "too big".
    2046             :  * Otherwise, is a constant times worse; for p = 3 (mod 4), is about 3 times worse,
    2047             :  * and in average is about 2 or 2.5 times worse. But try both algorithms for
    2048             :  * S(n) = (2^n+3)^2-8 with n = 750, 771, 779, 790, 874, 1176, 1728, 2604, etc.
    2049             :  *
    2050             :  * If X^2 := t^2 - a  is not a square in F_p (so X is in F_p^2), then
    2051             :  *   (t+X)^(p+1) = (t-X)(t+X) = a,   hence  sqrt(a) = (t+X)^((p+1)/2)  in F_p^2.
    2052             :  * If (a|p)=1, then sqrt(a) is in F_p.
    2053             :  * cf: LNCS 2286, pp 430-434 (2002)  [Gonzalo Tornaria] */
    2054             : 
    2055             : /* compute y^2, y = y[1] + y[2] X */
    2056             : static GEN
    2057         449 : sqrt_Cipolla_sqr(void *data, GEN y)
    2058             : {
    2059         449 :   GEN u = gel(y,1), v = gel(y,2), p = gel(data,2), n = gel(data,3);
    2060         449 :   GEN u2 = sqri(u), v2 = sqri(v);
    2061         449 :   v = subii(sqri(addii(v,u)), addii(u2,v2));
    2062         449 :   u = addii(u2, mulii(v2,n));
    2063             :   /* NOT mkvec2: must be gerepileupto-able */
    2064         449 :   retmkvec2(modii(u,p), modii(v,p));
    2065             : }
    2066             : /* compute (t+X) y^2 */
    2067             : static GEN
    2068          23 : sqrt_Cipolla_msqr(void *data, GEN y)
    2069             : {
    2070          23 :   GEN u = gel(y,1), v = gel(y,2), a = gel(data,1), p = gel(data,2), gt = gel(data,4);
    2071          23 :   ulong t = gt[2];
    2072          23 :   GEN d = addii(u, mului(t,v)), d2= sqri(d);
    2073          23 :   GEN b = remii(mulii(a,v), p);
    2074          23 :   u = subii(mului(t,d2), mulii(b,addii(u,d)));
    2075          23 :   v = subii(d2, mulii(b,v));
    2076             :   /* NOT mkvec2: must be gerepileupto-able */
    2077          23 :   retmkvec2(modii(u,p), modii(v,p));
    2078             : }
    2079             : /* assume a reduced mod p [ otherwise correct but inefficient ] */
    2080             : static GEN
    2081           8 : sqrt_Cipolla(GEN a, GEN p)
    2082             : {
    2083             :   pari_sp av1;
    2084             :   GEN u, v, n, y, pov2;
    2085             :   ulong t;
    2086             : 
    2087           8 :   if (kronecker(a, p) < 0) return NULL;
    2088           8 :   pov2 = shifti(p,-1);
    2089           8 :   if (cmpii(a,pov2) > 0) a = subii(a,p); /* center: avoid multiplying by huge base*/
    2090             : 
    2091           8 :   av1 = avma;
    2092          41 :   for(t=1; ; t++)
    2093             :   {
    2094          41 :     n = subsi((long)(t*t), a);
    2095          41 :     if (kronecker(n, p) < 0) break;
    2096          33 :     avma = av1;
    2097          33 :   }
    2098             : 
    2099             :   /* compute (t+X)^((p-1)/2) =: u+vX */
    2100           8 :   u = utoipos(t);
    2101           8 :   y = gen_pow_fold(mkvec2(u, gen_1), pov2, mkvec4(a,p,n,u),
    2102             :                          sqrt_Cipolla_sqr, sqrt_Cipolla_msqr);
    2103             :   /* Now u+vX = (t+X)^((p-1)/2); thus
    2104             :    *   (u+vX)(t+X) = sqrt(a) + 0 X
    2105             :    * Whence,
    2106             :    *   sqrt(a) = (u+vt)t - v*a
    2107             :    *   0       = (u+vt)
    2108             :    * Thus a square root is v*a */
    2109             : 
    2110           8 :   v = Fp_mul(gel(y, 2), a, p);
    2111           8 :   if (cmpii(v,pov2) > 0) v = subii(p,v);
    2112           8 :   return v;
    2113             : }
    2114             : 
    2115             : #define sqrmod(x,p) (remii(sqri(x),p))
    2116             : 
    2117             : /* Tonelli-Shanks. Assume p is prime and return NULL if (a,p) = -1. */
    2118             : GEN
    2119     2336851 : Fp_sqrt(GEN a, GEN p)
    2120             : {
    2121     2336851 :   pari_sp av = avma, av1;
    2122             :   long i, k, e;
    2123             :   GEN p1, q, v, y, w, m;
    2124             : 
    2125     2336851 :   if (typ(a) != t_INT) pari_err_TYPE("Fp_sqrt",a);
    2126     2336851 :   if (typ(p) != t_INT) pari_err_TYPE("Fp_sqrt",p);
    2127     2336851 :   if (signe(p) <= 0 || equali1(p)) pari_err_PRIME("Fp_sqrt",p);
    2128     2336851 :   if (lgefint(p) == 3)
    2129             :   {
    2130     2326741 :     ulong pp = uel(p,2), u = Fl_sqrt(umodiu(a, pp), pp);
    2131     2326727 :     if (u == ~0UL) return NULL;
    2132     2326685 :     return utoi(u);
    2133             :   }
    2134             : 
    2135       10110 :   p1 = addsi(-1,p); e = vali(p1);
    2136       10110 :   a = modii(a, p);
    2137             : 
    2138             :   /* On average, the algorithm of Cipolla is better than the algorithm of
    2139             :    * Tonelli and Shanks if and only if e(e-1)>8*log2(n)+20
    2140             :    * see LNCS 2286 pp 430 [GTL] */
    2141       10110 :   if (e*(e-1) > 20 + 8 * expi(p))
    2142             :   {
    2143           8 :     v = sqrt_Cipolla(a,p);
    2144           8 :     if (!v) { avma = av; return NULL; }
    2145           8 :     return gerepileuptoint(av,v);
    2146             :   }
    2147             : 
    2148       10102 :   if (e == 0) /* p = 2 */
    2149             :   {
    2150           0 :     avma = av;
    2151           0 :     if (!equaliu(p,2)) pari_err_PRIME("Fp_sqrt [modulus]",p);
    2152           0 :     if (!signe(a) || !mod2(a)) return gen_0;
    2153           0 :     return gen_1;
    2154             :   }
    2155       10102 :   q = shifti(p1,-e); /* q = (p-1)/2^oo is odd */
    2156       10102 :   if (e == 1) y = p1;
    2157             :   else /* look for an odd power of a primitive root */
    2158        8326 :     for (k=2; ; k++)
    2159             :     { /* loop terminates for k < p (even if p composite) */
    2160             : 
    2161        8326 :       i = krosi(k,p);
    2162        8326 :       if (i >= 0)
    2163             :       {
    2164        3131 :         if (i) continue;
    2165           0 :         pari_err_PRIME("Fp_sqrt [modulus]",p);
    2166             :       }
    2167        5195 :       av1 = avma;
    2168        5195 :       y = m = Fp_pow(utoipos((ulong)k),q,p);
    2169       12826 :       for (i=1; i<e; i++)
    2170        7631 :         if (gequal1(m = sqrmod(m,p))) break;
    2171        5195 :       if (i == e) break; /* success */
    2172           0 :       avma = av1;
    2173        3131 :     }
    2174             : 
    2175       10102 :   p1 = Fp_pow(a, shifti(q,-1), p); /* a ^ [(q-1)/2] */
    2176       10102 :   if (!signe(p1)) { avma=av; return gen_0; }
    2177       10095 :   v = Fp_mul(a, p1, p);
    2178       10095 :   w = Fp_mul(v, p1, p);
    2179       23990 :   while (!equali1(w))
    2180             :   { /* a*w = v^2, y primitive 2^e-th root of 1
    2181             :        a square --> w even power of y, hence w^(2^(e-1)) = 1 */
    2182        3810 :     p1 = sqrmod(w,p);
    2183        3810 :     for (k=1; !equali1(p1) && k < e; k++) p1 = sqrmod(p1,p);
    2184        3810 :     if (k == e) { avma=av; return NULL; } /* p composite or (a/p) != 1 */
    2185             :     /* w ^ (2^k) = 1 --> w = y ^ (u * 2^(e-k)), u odd */
    2186        3800 :     p1 = y;
    2187        3800 :     for (i=1; i < e-k; i++) p1 = sqrmod(p1,p);
    2188        3800 :     y = sqrmod(p1, p); e = k;
    2189        3800 :     w = Fp_mul(y, w, p);
    2190        3800 :     v = Fp_mul(v, p1, p);
    2191        3800 :     if (gc_needed(av,1))
    2192             :     {
    2193           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"Fp_sqrt");
    2194           0 :       gerepileall(av,3, &y,&w,&v);
    2195             :     }
    2196             :   }
    2197       10085 :   av1 = avma;
    2198       10085 :   p1 = subii(p,v); if (cmpii(v,p1) > 0) v = p1; else avma = av1;
    2199       10085 :   return gerepileuptoint(av, v);
    2200             : }
    2201             : 
    2202             : /*********************************************************************/
    2203             : /**                                                                 **/
    2204             : /**                        GCD & BEZOUT                             **/
    2205             : /**                                                                 **/
    2206             : /*********************************************************************/
    2207             : 
    2208             : GEN
    2209     3385341 : lcmii(GEN x, GEN y)
    2210             : {
    2211             :   pari_sp av;
    2212             :   GEN a, b;
    2213     3385341 :   if (!signe(x) || !signe(y)) return gen_0;
    2214     3385341 :   av = avma;
    2215     3385341 :   a = gcdii(x,y); if (!equali1(a)) y = diviiexact(y,a);
    2216     3385341 :   b = mulii(x,y); setabssign(b); return gerepileuptoint(av, b);
    2217             : }
    2218             : 
    2219             : /*********************************************************************/
    2220             : /**                                                                 **/
    2221             : /**                      CHINESE REMAINDERS                         **/
    2222             : /**                                                                 **/
    2223             : /*********************************************************************/
    2224             : 
    2225             : /* Chinese Remainder Theorem.  x and y must have the same type (integermod,
    2226             :  * polymod, or polynomial/vector/matrix recursively constructed with these
    2227             :  * as coefficients). Creates (with the same type) a z in the same residue
    2228             :  * class as x and the same residue class as y, if it is possible.
    2229             :  *
    2230             :  * We also allow (during recursion) two identical objects even if they are
    2231             :  * not integermod or polymod. For example:
    2232             :  *
    2233             :  * ? x = [1, Mod(5, 11), Mod(X + Mod(2, 7), X^2 + 1)];
    2234             :  * ? y = [1, Mod(7, 17), Mod(X + Mod(0, 3), X^2 + 1)];
    2235             :  * ? chinese(x, y)
    2236             :  * %3 = [1, Mod(16, 187), Mod(X + mod(9, 21), X^2 + 1)] */
    2237             : 
    2238             : static GEN
    2239      101761 : gen_chinese(GEN x, GEN(*f)(GEN,GEN))
    2240             : {
    2241      101761 :   GEN z = gassoc_proto(f,x,NULL);
    2242      101754 :   if (z == gen_1) retmkintmod(gen_0,gen_1);
    2243      101740 :   return z;
    2244             : }
    2245             : 
    2246             : /* x t_INTMOD, y t_POLMOD; promote x to t_POLMOD mod Pol(x.mod) then
    2247             :  * call chinese: makes Mod(0,1) a better "neutral" element */
    2248             : static GEN
    2249          21 : chinese_intpol(GEN x,GEN y)
    2250             : {
    2251          21 :   pari_sp av = avma;
    2252          21 :   GEN z = mkpolmod(gel(x,2), scalarpol_shallow(gel(x,1), varn(gel(y,1))));
    2253          21 :   return gerepileupto(av, chinese(z, y));
    2254             : }
    2255             : 
    2256             : GEN
    2257          49 : chinese1(GEN x) { return gen_chinese(x,chinese); }
    2258             : 
    2259             : GEN
    2260       16520 : chinese(GEN x, GEN y)
    2261             : {
    2262             :   pari_sp av;
    2263       16520 :   long tx = typ(x), ty;
    2264             :   GEN z,p1,p2,d,u,v;
    2265             : 
    2266       16520 :   if (!y) return chinese1(x);
    2267       16471 :   if (gidentical(x,y)) return gcopy(x);
    2268       16464 :   ty = typ(y);
    2269       16464 :   if (tx == ty) switch(tx)
    2270             :   {
    2271             :     case t_POLMOD:
    2272             :     {
    2273          28 :       GEN A = gel(x,1), B = gel(y,1);
    2274          28 :       GEN a = gel(x,2), b = gel(y,2);
    2275          28 :       if (varn(A)!=varn(B)) pari_err_VAR("chinese",A,B);
    2276          28 :       if (RgX_equal(A,B)) retmkpolmod(chinese(a,b), gcopy(A)); /*same modulus*/
    2277          28 :       av = avma;
    2278          28 :       d = RgX_extgcd(A,B,&u,&v);
    2279          28 :       p2 = gsub(b, a);
    2280          28 :       if (!gequal0(gmod(p2, d))) break;
    2281          28 :       p1 = gdiv(A,d);
    2282          28 :       p2 = gadd(a, gmul(gmul(u,p1), p2));
    2283             : 
    2284          28 :       z = cgetg(3, t_POLMOD);
    2285          28 :       gel(z,1) = gmul(p1,B);
    2286          28 :       gel(z,2) = gmod(p2,gel(z,1));
    2287          28 :       return gerepileupto(av, z);
    2288             :     }
    2289             :     case t_INTMOD:
    2290             :     {
    2291       16401 :       GEN A = gel(x,1), B = gel(y,1);
    2292       16401 :       GEN a = gel(x,2), b = gel(y,2), c, d, C, U;
    2293       16401 :       z = cgetg(3,t_INTMOD);
    2294       16401 :       Z_chinese_pre(A, B, &C, &U, &d);
    2295       16401 :       c = Z_chinese_post(a, b, C, U, d);
    2296       16401 :       if (!c) pari_err_OP("chinese", x,y);
    2297       16401 :       gel(z,1) = icopy_avma(C, (pari_sp)z);
    2298       16401 :       gel(z,2) = icopy_avma(c, (pari_sp)gel(z,1));
    2299       16401 :       avma = (pari_sp)gel(z,2); return z;
    2300             :     }
    2301             :     case t_POL:
    2302             :     {
    2303           7 :       long i, lx = lg(x), ly = lg(y);
    2304           7 :       if (varn(x) != varn(y)) break;
    2305           7 :       if (lx < ly) { swap(x,y); lswap(lx,ly); }
    2306           7 :       z = cgetg(lx, t_POL); z[1] = x[1];
    2307           7 :       for (i=2; i<ly; i++) gel(z,i) = chinese(gel(x,i),gel(y,i));
    2308           7 :       for (   ; i<lx; i++) gel(z,i) = gcopy(gel(x,i));
    2309           7 :       return z;
    2310             :     }
    2311             : 
    2312             :     case t_VEC: case t_COL: case t_MAT:
    2313             :     {
    2314             :       long i, lx;
    2315           7 :       z = cgetg_copy(x, &lx); if (lx!=lg(y)) break;
    2316           7 :       for (i=1; i<lx; i++) gel(z,i) = chinese(gel(x,i),gel(y,i));
    2317           7 :       return z;
    2318             :     }
    2319             :   }
    2320          21 :   if (tx == t_POLMOD && ty == t_INTMOD) return chinese_intpol(y,x);
    2321           7 :   if (ty == t_POLMOD && tx == t_INTMOD) return chinese_intpol(x,y);
    2322           0 :   pari_err_OP("chinese",x,y);
    2323           0 :   return NULL; /* not reached */
    2324             : }
    2325             : 
    2326             : /* init chinese(Mod(.,A), Mod(.,B)) */
    2327             : void
    2328      433741 : Z_chinese_pre(GEN A, GEN B, GEN *pC, GEN *pU, GEN *pd)
    2329             : {
    2330      433741 :   GEN u, d = bezout(A,B,&u,NULL); /* U = u(A/d), u(A/d) + v(B/d) = 1 */
    2331      433741 :   GEN t = diviiexact(A,d);
    2332      433741 :   *pU = mulii(u, t);
    2333      433741 :   *pC = mulii(t, B);
    2334      433741 :   if (pd) *pd = d;
    2335      433741 : }
    2336             : /* Assume C = lcm(A, B), U = 0 mod (A/d), U = 1 mod (B/d), a = b mod d,
    2337             :  * where d = gcd(A,B) or NULL, return x = a (mod A), b (mod B).
    2338             :  * If d not NULL, check wether a = b mod d. */
    2339             : GEN
    2340      590553 : Z_chinese_post(GEN a, GEN b, GEN C, GEN U, GEN d)
    2341             : {
    2342             :   GEN b_a;
    2343      590553 :   if (!signe(a))
    2344             :   {
    2345      437544 :     if (d && remii(b, d) != gen_0) return NULL;
    2346      437544 :     return Fp_mul(b, U, C);
    2347             :   }
    2348      153009 :   b_a = subii(b,a);
    2349      153009 :   if (d && remii(b_a, d) != gen_0) return NULL;
    2350      153009 :   return modii(addii(a, mulii(U, b_a)), C);
    2351             : }
    2352             : GEN
    2353        2142 : Z_chinese(GEN a, GEN b, GEN A, GEN B)
    2354             : {
    2355        2142 :   pari_sp av = avma;
    2356        2142 :   GEN C, U; Z_chinese_pre(A, B, &C, &U, NULL);
    2357        2142 :   return gerepileuptoint(av, Z_chinese_post(a,b, C, U, NULL));
    2358             : }
    2359             : GEN
    2360      415142 : Z_chinese_all(GEN a, GEN b, GEN A, GEN B, GEN *pC)
    2361             : {
    2362      415142 :   GEN U; Z_chinese_pre(A, B, pC, &U, NULL);
    2363      415142 :   return Z_chinese_post(a,b, *pC, U, NULL);
    2364             : }
    2365             : 
    2366             : /* return lift(chinese(a mod A, b mod B))
    2367             :  * assume(A,B)=1, a,b,A,B integers and C = A*B */
    2368             : GEN
    2369       83153 : Z_chinese_coprime(GEN a, GEN b, GEN A, GEN B, GEN C)
    2370             : {
    2371       83153 :   pari_sp av = avma;
    2372       83153 :   GEN U = mulii(Fp_inv(A,B), A);
    2373       83153 :   return gerepileuptoint(av, Z_chinese_post(a,b,C,U, NULL));
    2374             : }
    2375             : 
    2376             : /* chinese1 for coprime moduli in Z */
    2377             : static GEN
    2378       73393 : chinese1_coprime_Z_aux(GEN x, GEN y)
    2379             : {
    2380       73393 :   GEN z = cgetg(3, t_INTMOD);
    2381       73393 :   GEN A = gel(x,1), a = gel(x, 2);
    2382       73393 :   GEN B = gel(y,1), b = gel(y, 2), C = mulii(A,B);
    2383       73393 :   pari_sp av = avma;
    2384       73393 :   GEN U = mulii(Fp_inv(A,B), A);
    2385       73393 :   gel(z,2) = gerepileuptoint(av, Z_chinese_post(a,b,C,U, NULL));
    2386       73393 :   gel(z,1) = C; return z;
    2387             : }
    2388             : GEN
    2389      101712 : chinese1_coprime_Z(GEN x) {return gen_chinese(x,chinese1_coprime_Z_aux);}
    2390             : 
    2391             : /*********************************************************************/
    2392             : /**                                                                 **/
    2393             : /**                    MODULAR EXPONENTIATION                       **/
    2394             : /**                                                                 **/
    2395             : /*********************************************************************/
    2396             : 
    2397             : /* xa, ya = t_VECSMALL */
    2398             : GEN
    2399      155025 : ZV_producttree(GEN xa)
    2400             : {
    2401      155025 :   long n = lg(xa)-1;
    2402      155025 :   long m = n==1 ? 1: expu(n-1)+1;
    2403      155024 :   GEN T = cgetg(m+1, t_VEC), t;
    2404             :   long i, j, k;
    2405      155022 :   t = cgetg(((n+1)>>1)+1, t_VEC);
    2406      155023 :   if (typ(xa)==t_VECSMALL)
    2407             :   {
    2408      186773 :     for (j=1, k=1; k<n; j++, k+=2)
    2409      135133 :       gel(t, j) = muluu(xa[k], xa[k+1]);
    2410       51640 :     if (k==n) gel(t, j) = utoi(xa[k]);
    2411             :   } else {
    2412      279161 :     for (j=1, k=1; k<n; j++, k+=2)
    2413      175778 :       gel(t, j) = mulii(gel(xa,k), gel(xa,k+1));
    2414      103383 :     if (k==n) gel(t, j) = icopy(gel(xa,k));
    2415             :   }
    2416      155023 :   gel(T,1) = t;
    2417      279926 :   for (i=2; i<=m; i++)
    2418             :   {
    2419      124903 :     GEN u = gel(T, i-1);
    2420      124903 :     long n = lg(u)-1;
    2421      124903 :     t = cgetg(((n+1)>>1)+1, t_VEC);
    2422      331552 :     for (j=1, k=1; k<n; j++, k+=2)
    2423      206649 :       gel(t, j) = mulii(gel(u, k), gel(u, k+1));
    2424      124903 :     if (k==n) gel(t, j) = gel(u, k);
    2425      124903 :     gel(T, i) = t;
    2426             :   }
    2427      155023 :   return T;
    2428             : }
    2429             : 
    2430             : static GEN
    2431      844686 : Z_ZV_mod_tree(GEN P, GEN xa, GEN T)
    2432             : {
    2433             :   long i,j,k;
    2434      844686 :   long m = lg(T)-1, n = lg(xa)-1;
    2435             :   GEN t;
    2436      844686 :   GEN Tp = cgetg(m+1, t_VEC);
    2437      844667 :   gel(Tp, m) = mkvec(P);
    2438     1502256 :   for (i=m-1; i>=1; i--)
    2439             :   {
    2440      657608 :     GEN u = gel(T, i);
    2441      657608 :     GEN v = gel(Tp, i+1);
    2442      657608 :     long n = lg(u)-1;
    2443      657608 :     t = cgetg(n+1, t_VEC);
    2444     1611618 :     for (j=1, k=1; k<n; j++, k+=2)
    2445             :     {
    2446      953982 :       gel(t, k)   = modii(gel(v, j), gel(u, k));
    2447      953998 :       gel(t, k+1) = modii(gel(v, j), gel(u, k+1));
    2448             :     }
    2449      657636 :     if (k==n) gel(t, k) = gel(v, j);
    2450      657636 :     gel(Tp, i) = t;
    2451             :   }
    2452             :   {
    2453      844648 :     GEN u = gel(T, i+1);
    2454      844648 :     GEN v = gel(Tp, i+1);
    2455      844648 :     long l = lg(u)-1;
    2456      844648 :     if (typ(xa)==t_VECSMALL)
    2457             :     {
    2458      689622 :       GEN R = cgetg(n+1, t_VECSMALL);
    2459     2126632 :       for (j=1, k=1; j<=l; j++, k+=2)
    2460             :       {
    2461     1436957 :         uel(R,k) = umodiu(gel(v, j), xa[k]);
    2462     1436994 :         if (k < n)
    2463     1185942 :           uel(R,k+1) = umodiu(gel(v, j), xa[k+1]);
    2464             :       }
    2465      689675 :       return R;
    2466             :     }
    2467             :     else
    2468             :     {
    2469      155026 :       GEN R = cgetg(n+1, t_VEC);
    2470      516705 :       for (j=1, k=1; j<=l; j++, k+=2)
    2471             :       {
    2472      361675 :         gel(R,k) = modii(gel(v, j), gel(xa,k));
    2473      361677 :         if (k < n)
    2474      310914 :           gel(R,k+1) = modii(gel(v, j), gel(xa,k+1));
    2475             :       }
    2476      155030 :       return R;
    2477             :     }
    2478             :   }
    2479             : }
    2480             : 
    2481             : static GEN
    2482     2261122 : ZV_polint_tree(GEN T, GEN R, GEN xa, GEN ya)
    2483             : {
    2484     2261122 :   long m = lg(T)-1, n = lg(ya)-1;
    2485             :   long i,j,k;
    2486     2261122 :   GEN Tp = cgetg(m+1, t_VEC);
    2487     2250854 :   GEN M = gel(T, 1);
    2488     2250854 :   GEN t = cgetg(lg(M), t_VEC);
    2489     2297719 :   if (typ(xa)==t_VECSMALL)
    2490             :   {
    2491    18015834 :     for (j=1, k=1; k<n; j++, k+=2)
    2492             :     {
    2493    15858914 :       pari_sp av = avma;
    2494    15858914 :       GEN a = mului(ya[k], gel(R,k)), b = mului(ya[k+1], gel(R,k+1));
    2495    15693833 :       GEN tj = modii(addii(mului(xa[k],b), mului(xa[k+1],a)), gel(M,j));
    2496    15680550 :       gel(t, j) = gerepileuptoint(av, tj);
    2497             :     }
    2498     2156920 :     if (k==n) gel(t, j) = modii(mului(ya[k], gel(R,k)), gel(M, j));
    2499             :   } else
    2500             :   {
    2501      279161 :     for (j=1, k=1; k<n; j++, k+=2)
    2502             :     {
    2503      175778 :       pari_sp av = avma;
    2504      175778 :       GEN a = mulii(gel(ya,k), gel(R,k)), b = mulii(gel(ya,k+1), gel(R,k+1));
    2505      175778 :       GEN tj = modii(addii(mulii(gel(xa,k),b), mulii(gel(xa,k+1),a)), gel(M,j));
    2506      175778 :       gel(t, j) = gerepileuptoint(av, tj);
    2507             :     }
    2508      103383 :     if (k==n) gel(t, j) = modii(mulii(gel(ya,k), gel(R,k)), gel(M, j));
    2509             :   }
    2510     2259822 :   gel(Tp, 1) = t;
    2511     8449333 :   for (i=2; i<=m; i++)
    2512             :   {
    2513     6200357 :     GEN u = gel(T, i-1), M = gel(T, i);
    2514     6200357 :     GEN t = cgetg(lg(M), t_VEC);
    2515     6232009 :     GEN v = gel(Tp, i-1);
    2516     6232009 :     long n = lg(v)-1;
    2517    20576399 :     for (j=1, k=1; k<n; j++, k+=2)
    2518             :     {
    2519    14386888 :       pari_sp av = avma;
    2520    57547552 :       gel(t, j) = gerepileuptoint(av, modii(addii(mulii(gel(u, k), gel(v, k+1)),
    2521    43160664 :             mulii(gel(u, k+1), gel(v, k))), gel(M, j)));
    2522             :     }
    2523     6189511 :     if (k==n) gel(t, j) = gel(v, k);
    2524     6189511 :     gel(Tp, i) = t;
    2525             :   }
    2526     2248976 :   return gmael(Tp,m,1);
    2527             : }
    2528             : 
    2529             : static GEN
    2530     2118559 : ZV_polint_center_tree(GEN T, GEN R, GEN xa, GEN ya, GEN m2)
    2531             : {
    2532     2118559 :   GEN mod = gmael(T, lg(T)-1, 1);
    2533     2118559 :   GEN a = ZV_polint_tree(T, R, xa, ya);
    2534     2096622 :   return Fp_center(a, mod, m2);
    2535             : }
    2536             : 
    2537             : static GEN
    2538       40819 : ncV_polint_center_tree(GEN T, GEN R, GEN xa, GEN Va, GEN m2)
    2539             : {
    2540       40819 :   long i, j, l = lg(gel(Va,1)), n = lg(xa);
    2541       40819 :   GEN V = cgetg(l, t_COL);
    2542     2162301 :   for(i=1; i < l; i++)
    2543             :   {
    2544     2121296 :     pari_sp av = avma;
    2545     2121296 :     GEN ya = cgetg(n, t_VECSMALL);
    2546    34888557 :     for(j=1; j < n; j++)
    2547    32755784 :       ya[j] = mael(Va,j,i);
    2548     2132773 :     gel(V,i) = gerepilecopy(av, ZV_polint_center_tree(T, R, xa, ya, m2));
    2549             :   }
    2550       41005 :   return V;
    2551             : }
    2552             : 
    2553             : GEN
    2554       38030 : nmV_polint_center_tree_worker(GEN Va, GEN T, GEN R, GEN xa, GEN m2)
    2555             : {
    2556       38030 :   return ncV_polint_center_tree(T, R, xa, Va, m2);
    2557             : }
    2558             : 
    2559             : static GEN
    2560        2090 : nmV_polint_center_tree(GEN T, GEN R, GEN xa, GEN Ma, GEN m2)
    2561             : {
    2562        2090 :   long i, j, l = lg(gel(Ma,1)), n = lg(xa);
    2563        2090 :   long pending = 0, workid, cnt = 0;
    2564             :   struct pari_mt pt;
    2565             :   GEN worker, done, va, M;
    2566        2090 :   GEN ya = cgetg(n, t_VEC);
    2567        2090 :   worker = snm_closure(is_entry("_polint_worker"), mkvec4(T, R, xa, m2));
    2568        2090 :   va = mkvec(gen_0);
    2569        2090 :   M = cgetg(l, t_MAT);
    2570        2090 :   if (DEBUGLEVEL>2) err_printf("Start parallel Chinese remainder: ");
    2571        2090 :   mt_queue_start(&pt, worker);
    2572       43765 :   for (i=1; i<l || pending; i++)
    2573             :   {
    2574      611825 :     for(j=1; j < n; j++)
    2575      570150 :       gel(ya,j) = gmael(Ma,j,i);
    2576       41675 :     gel(va, 1) = ya;
    2577       41675 :     mt_queue_submit(&pt, i, i<l? va: NULL);
    2578       41675 :     done = mt_queue_get(&pt, &workid, &pending);
    2579       41675 :     if (done)
    2580             :     {
    2581       38256 :       gel(M,workid) = done;
    2582       38256 :       if (DEBUGLEVEL>2) err_printf("%ld%% ",(++cnt)*100/(l-1));
    2583             :     }
    2584             :   }
    2585        2090 :   if (DEBUGLEVEL>2) err_printf("\n");
    2586        2090 :   mt_queue_end(&pt);
    2587        2090 :   return M;
    2588             : }
    2589             : 
    2590             : GEN
    2591           0 : Z_ZV_mod(GEN P, GEN xa)
    2592             : {
    2593           0 :   pari_sp av = avma;
    2594           0 :   GEN T = ZV_producttree(xa);
    2595           0 :   return gerepileuptoleaf(av, Z_ZV_mod_tree(P, xa, T));
    2596             : }
    2597             : 
    2598             : GEN
    2599           0 : Z_nv_mod(GEN P, GEN xa)
    2600             : {
    2601           0 :   return Z_ZV_mod(P, xa);
    2602             : }
    2603             : 
    2604             : GEN
    2605       58813 : ZX_nv_mod_tree(GEN P, GEN xa, GEN T)
    2606             : {
    2607       58813 :   long i, j, l = lg(P), n = lg(xa)-1;
    2608       58813 :   GEN V = cgetg(n+1, t_VEC);
    2609      252228 :   for (j=1; j <= n; j++)
    2610             :   {
    2611      193400 :     gel(V, j) = cgetg(l, t_VECSMALL);
    2612      193399 :     mael(V, j, 1) = P[1]&VARNBITS;
    2613             :   }
    2614      748499 :   for (i=2; i < l; i++)
    2615             :   {
    2616      689681 :     GEN v = Z_ZV_mod_tree(gel(P, i), xa, T);
    2617     3312615 :     for (j=1; j <= n; j++)
    2618     2622944 :       mael(V, j, i) = v[j];
    2619             :   }
    2620      252227 :   for (j=1; j <= n; j++)
    2621      193411 :     (void) Flx_renormalize(gel(V, j), l);
    2622       58816 :   return V;
    2623             : }
    2624             : 
    2625             : static GEN
    2626      155028 : ZV_sqr(GEN z)
    2627             : {
    2628      155028 :   long i,l = lg(z);
    2629      155028 :   GEN x = cgetg(l, t_VEC);
    2630      155027 :   if (typ(z)==t_VECSMALL)
    2631       51644 :     for (i=1; i<l; i++) gel(x,i) = sqru(z[i]);
    2632             :   else
    2633      103383 :     for (i=1; i<l; i++) gel(x,i) = sqri(gel(z,i));
    2634      155028 :   return x;
    2635             : }
    2636             : 
    2637             : static GEN
    2638     1032724 : ZT_sqr(GEN z)
    2639             : {
    2640     1032724 :   if (typ(z) == t_INT)
    2641      597775 :     return sqri(z);
    2642             :   else
    2643             :   {
    2644      434949 :     long i,l = lg(z);
    2645      434949 :     GEN x = cgetg(l, t_VEC);
    2646      434950 :     for (i=1; i<l; i++) gel(x,i) = ZT_sqr(gel(z,i));
    2647      434958 :     return x;
    2648             :   }
    2649             : }
    2650             : 
    2651             : static GEN
    2652      155028 : ZV_invdivexact(GEN y, GEN x)
    2653             : {
    2654      155028 :   long i, l = lg(y);
    2655      155028 :   GEN z = cgetg(l,t_VEC);
    2656      155028 :   if (typ(x)==t_VECSMALL)
    2657      338715 :     for (i=1; i<l; i++)
    2658             :     {
    2659      287073 :       pari_sp av = avma;
    2660      287073 :       ulong a = Fl_inv(umodiu(diviuexact(gel(y,i),x[i]), x[i]), x[i]);
    2661      287079 :       avma = av;
    2662      287079 :       gel(z,i) = utoi(a);
    2663             :     }
    2664             :   else
    2665      488897 :     for (i=1; i<l; i++)
    2666      385514 :       gel(z,i) = Fp_inv(diviiexact(gel(y,i), gel(x,i)), gel(x,i));
    2667      155025 :   return z;
    2668             : }
    2669             : 
    2670             : static GEN
    2671      155023 : ZV_chinesetree(GEN T, GEN xa)
    2672             : {
    2673      155023 :   GEN T2 = ZT_sqr(T), xa2 = ZV_sqr(xa);
    2674      155028 :   GEN mod = gmael(T,lg(T)-1,1);
    2675      155028 :   return ZV_invdivexact(Z_ZV_mod_tree(mod, xa2, T2), xa);
    2676             : }
    2677             : 
    2678             : static GEN
    2679      155023 : gc_chinese(pari_sp av, GEN T, GEN a, GEN *pt_mod)
    2680             : {
    2681      155023 :   if (!pt_mod)
    2682        2090 :     return gerepileupto(av, a);
    2683             :   else
    2684             :   {
    2685      152933 :     GEN mod = gmael(T, lg(T)-1, 1);
    2686      152933 :     gerepileall(av, 2, &a, &mod);
    2687      152939 :     *pt_mod = mod;
    2688      152939 :     return a;
    2689             :   }
    2690             : }
    2691             : 
    2692             : GEN
    2693       46758 : ZV_chinese_tree(GEN A, GEN P, GEN T, GEN *pt_mod)
    2694             : {
    2695       46758 :   pari_sp av = avma;
    2696       46758 :   GEN R = ZV_chinesetree(T, P);
    2697       46756 :   GEN a = ZV_polint_tree(T, R, P, A);
    2698       46756 :   return gc_chinese(av, T, a, pt_mod);
    2699             : }
    2700             : 
    2701             : GEN
    2702      103383 : ZV_chinese(GEN A, GEN P, GEN *pt_mod)
    2703             : {
    2704      103383 :   pari_sp av = avma;
    2705      103383 :   GEN T = ZV_producttree(P);
    2706      103383 :   GEN R = ZV_chinesetree(T, P);
    2707      103383 :   GEN a = ZV_polint_tree(T, R, P, A);
    2708      103383 :   return gc_chinese(av, T, a, pt_mod);
    2709             : }
    2710             : 
    2711             : GEN
    2712        2794 : ncV_chinese_center(GEN A, GEN P, GEN *pt_mod)
    2713             : {
    2714        2794 :   pari_sp av = avma;
    2715        2794 :   GEN T = ZV_producttree(P);
    2716        2794 :   GEN R = ZV_chinesetree(T, P);
    2717        2794 :   GEN m2 = shifti(gmael(T, lg(T)-1, 1), -1);
    2718        2794 :   GEN a = ncV_polint_center_tree(T, R, P, A, m2);
    2719        2794 :   return gc_chinese(av, T, a, pt_mod);
    2720             : }
    2721             : 
    2722             : GEN
    2723        2090 : nmV_chinese_center(GEN A, GEN P, GEN *pt_mod)
    2724             : {
    2725        2090 :   pari_sp av = avma;
    2726        2090 :   GEN T = ZV_producttree(P);
    2727        2090 :   GEN R = ZV_chinesetree(T, P);
    2728        2090 :   GEN m2 = shifti(gmael(T, lg(T)-1, 1), -1);
    2729        2090 :   GEN a = nmV_polint_center_tree(T, R, P, A, m2);
    2730        2090 :   return gc_chinese(av, T, a, pt_mod);
    2731             : }
    2732             : 
    2733             : /**********************************************************************
    2734             :  **                                                                  **
    2735             :  **                    Powering  over (Z/NZ)^*, small N              **
    2736             :  **                                                                  **
    2737             :  **********************************************************************/
    2738             : 
    2739             : ulong
    2740   116631589 : Fl_powu_pre(ulong x, ulong n0, ulong p, ulong pi)
    2741             : {
    2742             :   ulong y, z, n;
    2743   116631589 :   if (n0 <= 1)
    2744             :   { /* frequent special cases */
    2745     9158117 :     if (n0 == 1) return x;
    2746     3476504 :     if (n0 == 0) return 1;
    2747             :   }
    2748   107473472 :   if (x <= 1) return x; /* 0 or 1 */
    2749   106455937 :   y = 1; z = x; n = n0;
    2750             :   for(;;)
    2751             :   {
    2752  3089890736 :     if (n&1) y = Fl_mul_pre(y,z,p,pi);
    2753  3089400136 :     n>>=1; if (!n) return y;
    2754  2982945417 :     z = Fl_sqr_pre(z,p,pi);
    2755  2983434799 :   }
    2756             : }
    2757             : 
    2758             : ulong
    2759    61034605 : Fl_powu(ulong x, ulong n0, ulong p)
    2760             : {
    2761             :   ulong y, z, n;
    2762    61034605 :   if (n0 <= 2)
    2763             :   { /* frequent special cases */
    2764    16114459 :     if (n0 == 2) return Fl_sqr(x,p);
    2765     2123582 :     if (n0 == 1) return x;
    2766        2638 :     if (n0 == 0) return 1;
    2767             :   }
    2768    44920146 :   if (x <= 1) return x; /* 0 or 1 */
    2769    44882994 :   if (!SMALL_ULONG(p))
    2770    40785308 :     return Fl_powu_pre(x, n0, p, get_Fl_red(p));
    2771     4097686 :   y = 1; z = x; n = n0;
    2772             :   for(;;)
    2773             :   {
    2774    53118150 :     if (n&1) y = Fl_mul(y,z,p);
    2775    53136902 :     n>>=1; if (!n) return y;
    2776    49032912 :     z = Fl_sqr(z,p);
    2777    49020464 :   }
    2778             : }
    2779             : 
    2780             : /* Reduce data dependency to maximize internal parallelism */
    2781             : GEN
    2782     8284826 : Fl_powers_pre(ulong x, long n, ulong p, ulong pi)
    2783             : {
    2784             :   long i, k;
    2785     8284826 :   GEN powers = cgetg(n + 2, t_VECSMALL);
    2786     8285440 :   powers[1] = 1; if (n == 0) return powers;
    2787     8285440 :   powers[2] = x;
    2788    43907420 :   for (i = 3, k=2; i <= n; i+=2, k++)
    2789             :   {
    2790    35618423 :     powers[i] = Fl_mul_pre(powers[k], powers[k], p, pi);
    2791    35621481 :     powers[i+1] = Fl_mul_pre(powers[k], powers[k+1], p, pi);
    2792             :   }
    2793     8288997 :   if (i==n+1)
    2794     7352118 :     powers[i] = Fl_mul_pre(powers[k], powers[k], p, pi);
    2795     8289006 :   return powers;
    2796             : }
    2797             : 
    2798             : GEN
    2799        1520 : Fl_powers(ulong x, long n, ulong p)
    2800             : {
    2801        1520 :   return Fl_powers_pre(x, n, p, get_Fl_red(p));
    2802             : }
    2803             : 
    2804             : /**********************************************************************
    2805             :  **                                                                  **
    2806             :  **                    Powering  over (Z/NZ)^*, large N              **
    2807             :  **                                                                  **
    2808             :  **********************************************************************/
    2809             : 
    2810             : static GEN
    2811      387210 : Fp_dblsqr(GEN x, GEN N)
    2812             : {
    2813      387210 :   GEN z = shifti(Fp_sqr(x, N), 1);
    2814      387210 :   return cmpii(z, N) >= 0? subii(z, N): z;
    2815             : }
    2816             : 
    2817             : typedef struct muldata {
    2818             :   GEN (*sqr)(void * E, GEN x);
    2819             :   GEN (*mul)(void * E, GEN x, GEN y);
    2820             :   GEN (*mul2)(void * E, GEN x);
    2821             : } muldata;
    2822             : 
    2823             : /* modified Barrett reduction with one fold */
    2824             : /* See Fast Modular Reduction, W. Hasenplaugh, G. Gaubatz, V. Gopal, ARITH 18 */
    2825             : 
    2826             : static GEN
    2827        5166 : Fp_invmBarrett(GEN p, long s)
    2828             : {
    2829        5166 :   GEN R, Q = dvmdii(int2n(3*s),p,&R);
    2830        5166 :   return mkvec2(Q,R);
    2831             : }
    2832             : 
    2833             : /* a <= (N-1)^2, 2^(2s-2) <= N < 2^(2s). Return 0 <= r < N such that
    2834             :  * a = r (mod N) */
    2835             : static GEN
    2836      316949 : Fp_rem_mBarrett(GEN a, GEN B, long s, GEN N)
    2837             : {
    2838      316949 :   pari_sp av = avma;
    2839      316949 :   GEN P = gel(B, 1), Q = gel(B, 2); /* 2^(3s) = P N + Q, 0 <= Q < N */
    2840      316949 :   long t = expi(P)+1; /* 2^(t-1) <= P < 2^t */
    2841      316949 :   GEN u = shifti(a, -3*s), v = remi2n(a, 3*s); /* a = 2^(3s)u + v */
    2842      316949 :   GEN A = addii(v, mulii(Q,u)); /* 0 <= A < 2^(3s+1) */
    2843      316949 :   GEN q = shifti(mulii(shifti(A, t-3*s), P), -t); /* A/N - 4 < q <= A/N */
    2844      316949 :   GEN r = subii(A, mulii(q, N));
    2845      316949 :   GEN sr= subii(r,N);     /* 0 <= r < 4*N */
    2846      316949 :   if (signe(sr)<0) return gerepileuptoint(av, r);
    2847      304699 :   r=sr; sr = subii(r,N);  /* 0 <= r < 3*N */
    2848      304699 :   if (signe(sr)<0) return gerepileuptoint(av, r);
    2849       17439 :   r=sr; sr = subii(r,N);  /* 0 <= r < 2*N */
    2850       17439 :   return gerepileuptoint(av, signe(sr)>=0 ? sr:r);
    2851             : }
    2852             : 
    2853             : /* Montgomery reduction */
    2854             : 
    2855             : INLINE ulong
    2856      239719 : init_montdata(GEN N) { return (ulong) -invmod2BIL(mod2BIL(N)); }
    2857             : 
    2858             : struct montred
    2859             : {
    2860             :   GEN N;
    2861             :   ulong inv;
    2862             : };
    2863             : 
    2864             : /* Montgomery reduction */
    2865             : static GEN
    2866    11852515 : _sqr_montred(void * E, GEN x)
    2867             : {
    2868    11852515 :   struct montred * D = (struct montred *) E;
    2869    11852515 :   return red_montgomery(sqri(x), D->N, D->inv);
    2870             : }
    2871             : 
    2872             : /* Montgomery reduction */
    2873             : static GEN
    2874     1238052 : _mul_montred(void * E, GEN x, GEN y)
    2875             : {
    2876     1238052 :   struct montred * D = (struct montred *) E;
    2877     1238052 :   return red_montgomery(mulii(x, y), D->N, D->inv);
    2878             : }
    2879             : 
    2880             : static GEN
    2881      813330 : _mul2_montred(void * E, GEN x)
    2882             : {
    2883      813330 :   struct montred * D = (struct montred *) E;
    2884      813330 :   GEN z = shifti(_sqr_montred(E, x), 1);
    2885      813329 :   long l = lgefint(D->N);
    2886      813329 :   while (lgefint(z) > l) z = subii(z, D->N);
    2887      813330 :   return z;
    2888             : }
    2889             : 
    2890             : static GEN
    2891     3630050 : _sqr_remii(void* N, GEN x)
    2892     3630050 : { return remii(sqri(x), (GEN) N); }
    2893             : 
    2894             : static GEN
    2895      240726 : _mul_remii(void* N, GEN x, GEN y)
    2896      240726 : { return remii(mulii(x, y), (GEN) N); }
    2897             : 
    2898             : static GEN
    2899      382690 : _mul2_remii(void* N, GEN x)
    2900      382690 : { return Fp_dblsqr(x, (GEN) N); }
    2901             : 
    2902             : struct redbarrett
    2903             : {
    2904             :   GEN iM, N;
    2905             :   long s;
    2906             : };
    2907             : 
    2908             : static GEN
    2909      308546 : _sqr_remiibar(void *E, GEN x)
    2910             : {
    2911      308546 :   struct redbarrett * D = (struct redbarrett *) E;
    2912      308546 :   return Fp_rem_mBarrett(sqri(x), D->iM, D->s, D->N);
    2913             : }
    2914             : 
    2915             : static GEN
    2916        8403 : _mul_remiibar(void *E, GEN x, GEN y)
    2917             : {
    2918        8403 :   struct redbarrett * D = (struct redbarrett *) E;
    2919        8403 :   return Fp_rem_mBarrett(mulii(x, y), D->iM, D->s, D->N);
    2920             : }
    2921             : 
    2922             : static GEN
    2923        4520 : _mul2_remiibar(void *E, GEN x)
    2924             : {
    2925        4520 :   struct redbarrett * D = (struct redbarrett *) E;
    2926        4520 :   return Fp_dblsqr(x, D->N);
    2927             : }
    2928             : 
    2929             : static long
    2930      312511 : Fp_select_red(GEN *y, ulong k, GEN N, long lN, muldata *D, void **pt_E)
    2931             : {
    2932      312511 :   if (lN >= Fp_POW_BARRETT_LIMIT && (k==0 || ((double)k)*expi(*y) > 2 + expi(N)))
    2933             :   {
    2934        5166 :     struct redbarrett * E = (struct redbarrett *) stack_malloc(sizeof(struct redbarrett));
    2935        5166 :     D->sqr = &_sqr_remiibar;
    2936        5166 :     D->mul = &_mul_remiibar;
    2937        5166 :     D->mul2 = &_mul2_remiibar;
    2938        5166 :     E->N = N;
    2939        5166 :     E->s = 1+(expi(N)>>1);
    2940        5166 :     E->iM = Fp_invmBarrett(N, E->s);
    2941        5166 :     *pt_E = (void*) E;
    2942        5166 :     return 0;
    2943             :   }
    2944      307345 :   else if (mod2(N) && lN < Fp_POW_REDC_LIMIT)
    2945             :   {
    2946      239719 :     struct montred * E = (struct montred *) stack_malloc(sizeof(struct montred));
    2947      239719 :     *y = remii(shifti(*y, bit_accuracy(lN)), N);
    2948      239719 :     D->sqr = &_sqr_montred;
    2949      239719 :     D->mul = &_mul_montred;
    2950      239719 :     D->mul2 = &_mul2_montred;
    2951      239719 :     E->N = N;
    2952      239719 :     E->inv = init_montdata(N);
    2953      239719 :     *pt_E = (void*) E;
    2954      239719 :     return 1;
    2955             :   }
    2956             :   else
    2957             :   {
    2958       67626 :     D->sqr = &_sqr_remii;
    2959       67626 :     D->mul = &_mul_remii;
    2960       67626 :     D->mul2 = &_mul2_remii;
    2961       67626 :     *pt_E = (void*) N;
    2962       67626 :     return 0;
    2963             :   }
    2964             : }
    2965             : 
    2966             : GEN
    2967      860183 : Fp_powu(GEN A, ulong k, GEN N)
    2968             : {
    2969      860183 :   long lN = lgefint(N), sA;
    2970             :   int base_is_2, use_montgomery;
    2971             :   muldata D;
    2972             :   void *E;
    2973             :   pari_sp av;
    2974             : 
    2975      860183 :   if (lN == 3) {
    2976      117865 :     ulong n = uel(N,2);
    2977      117865 :     return utoi( Fl_powu(umodiu(A, n), k, n) );
    2978             :   }
    2979      742318 :   if (k <= 2)
    2980             :   { /* frequent special cases */
    2981      522859 :     if (k == 2) return Fp_sqr(A,N);
    2982      119056 :     if (k == 1) return A;
    2983           0 :     if (k == 0) return gen_1;
    2984             :   }
    2985      219459 :   sA = signe(A)==-1 && odd(k);
    2986      219459 :   base_is_2 = 0;
    2987      219459 :   if (lgefint(A) == 3) switch(A[2])
    2988             :   {
    2989         492 :     case 1: return sA ? gen_m1 : gen_1;
    2990       36762 :     case 2:  base_is_2 = 1; break;
    2991             :   }
    2992             : 
    2993             :   /* TODO: Move this out of here and use for general modular computations */
    2994      218967 :   av = avma;
    2995      218967 :   use_montgomery = Fp_select_red(&A, k, N, lN, &D, &E);
    2996      218967 :   if (base_is_2)
    2997       36762 :     A = gen_powu_fold_i(A, k, E, D.sqr, D.mul2);
    2998             :   else
    2999      182205 :     A = gen_powu_i(A, k, E, D.sqr, D.mul);
    3000      218967 :   if (use_montgomery)
    3001             :   {
    3002      195663 :     A = red_montgomery(A, N, ((struct montred *) E)->inv);
    3003      195663 :     if (cmpii(A, N) >= 0) A = subii(A,N);
    3004      195663 :     if (sA) A = subii(N, A);
    3005             :   }
    3006      218967 :   return gerepileuptoint(av, A);
    3007             : }
    3008             : 
    3009             : GEN
    3010       21000 : Fp_pows(GEN A, long k, GEN N)
    3011             : {
    3012       21000 :   if (lgefint(N) == 3) {
    3013        6973 :     ulong n = N[2];
    3014        6973 :     ulong a = umodiu(A, n);
    3015        6973 :     if (k < 0) {
    3016          14 :       a = Fl_inv(a, n);
    3017          14 :       k = -k;
    3018             :     }
    3019        6973 :     return utoi( Fl_powu(a, (ulong)k, n) );
    3020             :   }
    3021       14027 :   if (k < 0) { A = Fp_inv(A, N); k = -k; };
    3022       14027 :   return Fp_powu(A, (ulong)k, N);
    3023             : }
    3024             : 
    3025             : /* A^K mod N */
    3026             : GEN
    3027     2497787 : Fp_pow(GEN A, GEN K, GEN N)
    3028             : {
    3029     2497787 :   pari_sp av = avma;
    3030     2497787 :   long t,s, lN = lgefint(N), sA;
    3031             :   int base_is_2, use_montgomery;
    3032             :   GEN y;
    3033             :   muldata D;
    3034             :   void *E;
    3035             : 
    3036     2497787 :   s = signe(K);
    3037     2497787 :   if (!s)
    3038             :   {
    3039       14980 :     t = signe(remii(A,N)); avma = av;
    3040       14980 :     return t? gen_1: gen_0;
    3041             :   }
    3042     2482807 :   if (lN == 3)
    3043             :   {
    3044     2243076 :     ulong k, n = N[2], a = umodiu(A, n);
    3045     2243076 :     if (s < 0) a = Fl_inv(a, n);
    3046     2243020 :     if (a <= 1) return utoi(a); /* 0 or 1 */
    3047     2132632 :     if (lgefint(K) > 3)
    3048             :     { /* silly case : huge exponent, small modulus */
    3049          21 :       pari_warn(warner, "Mod(a,b)^n with n >> b : wasteful");
    3050          21 :       if (s > 0)
    3051             :       {
    3052          14 :         ulong d = ugcd(a, n);
    3053          14 :         if (d != 1)
    3054             :         { /* write n = n1 n2, with n2 maximal such that (n1,a) = 1 */
    3055           7 :           ulong n1 = ucoprime_part(n, d), n2 = n/n1;
    3056             : 
    3057           7 :           k = umodiu(K, eulerphiu(n1));
    3058             :           /* CRT: = a^K (mod n1), = 0 (mod n2)*/
    3059           7 :           return utoi( Fl_mul(Fl_powu(a, k, n1), n2 * Fl_inv(n2,n1), n) );
    3060             :         }
    3061             :         /* gcd(a,n) = 1 */
    3062           7 :         k = umodiu(K, eulerphiu(n));
    3063             :       }
    3064             :       else
    3065           7 :         k = umodiu(negi(K), eulerphiu(n));
    3066             :     }
    3067             :     else
    3068     2132611 :       k = uel(K,2);
    3069     2132625 :     return utoi(Fl_powu(a, k, n));
    3070             :   }
    3071             : 
    3072      239731 :   if (s < 0) y = Fp_inv(A,N);
    3073             :   else
    3074             :   {
    3075      239380 :     y = modii(A,N);
    3076      239380 :     if (!signe(y)) { avma = av; return gen_0; }
    3077             :   }
    3078      239724 :   if (lgefint(K) == 3) return gerepileuptoint(av, Fp_powu(y, K[2], N));
    3079             : 
    3080       93605 :   base_is_2 = 0;
    3081       93605 :   sA = signe(y)==-1 && mod2(K);
    3082       93605 :   if (lgefint(y) == 3) switch(y[2])
    3083             :   {
    3084          61 :     case 1: return sA ? gen_m1 : gen_1;
    3085       64728 :     case 2:  base_is_2 = 1; break;
    3086             :   }
    3087             : 
    3088             :   /* TODO: Move this out of here and use for general modular computations */
    3089       93544 :   use_montgomery = Fp_select_red(&y, 0UL, N, lN, &D, &E);
    3090       93544 :   if (base_is_2)
    3091       64728 :     y = gen_pow_fold_i(y, K, E, D.sqr, D.mul2);
    3092             :   else
    3093       28816 :     y = gen_pow_i(y, K, E, D.sqr, D.mul);
    3094       93543 :   if (use_montgomery)
    3095             :   {
    3096       44056 :     y = red_montgomery(y, N, ((struct montred *) E)->inv);
    3097       44056 :     if (cmpii(y,N) >= 0) y = subii(y,N);
    3098       44055 :     if (sA) y = subii(N, y);
    3099             :   }
    3100       93542 :   return gerepileuptoint(av,y);
    3101             : }
    3102             : 
    3103             : static GEN
    3104      166047 : _Fp_mul(void *E, GEN x, GEN y) { return Fp_mul(x,y,(GEN)E); }
    3105             : 
    3106             : static GEN
    3107         948 : _Fp_sqr(void *E, GEN x) { return Fp_sqr(x,(GEN)E); }
    3108             : 
    3109             : static GEN
    3110         104 : _Fp_one(void *E) { (void) E; return gen_1; }
    3111             : 
    3112             : GEN
    3113        1624 : Fp_powers(GEN x, long n, GEN p)
    3114             : {
    3115        1624 :   if (lgefint(p) == 3)
    3116        1520 :     return Flv_to_ZV(Fl_powers(umodiu(x, uel(p, 2)), n, uel(p, 2)));
    3117         104 :   return gen_powers(x, n, 1, (void*)p, _Fp_sqr, _Fp_mul, _Fp_one);
    3118             : }
    3119             : 
    3120             : static GEN
    3121      625233 : _Fp_pow(void *E, GEN x, GEN n) { return Fp_pow(x,n,(GEN)E); }
    3122             : 
    3123             : static GEN
    3124         413 : _Fp_rand(void *E) { return addis(randomi(subis((GEN)E,1)),1); }
    3125             : 
    3126             : static GEN Fp_easylog(void *E, GEN a, GEN g, GEN ord);
    3127             : 
    3128             : static const struct bb_group Fp_star={_Fp_mul,_Fp_pow,_Fp_rand,hash_GEN,
    3129             :                                       equalii,equali1,Fp_easylog};
    3130             : 
    3131             : static GEN
    3132      935783 : _Fp_red(void *E, GEN x) { return Fp_red(x, (GEN)E); }
    3133             : 
    3134             : static GEN
    3135      944352 : _Fp_add(void *E, GEN x, GEN y) { (void) E; return addii(x,y); }
    3136             : 
    3137             : static GEN
    3138      103113 : _Fp_neg(void *E, GEN x) { (void) E; return negi(x); }
    3139             : 
    3140             : static GEN
    3141     1181019 : _Fp_rmul(void *E, GEN x, GEN y) { (void) E; return mulii(x,y); }
    3142             : 
    3143             : static GEN
    3144       18323 : _Fp_inv(void *E, GEN x) { return Fp_inv(x,(GEN)E); }
    3145             : 
    3146             : static int
    3147      526253 : _Fp_equal0(GEN x) { return signe(x)==0; }
    3148             : 
    3149             : static GEN
    3150       82540 : _Fp_s(void *E, long x) { (void) E; return stoi(x); }
    3151             : 
    3152             : static const struct bb_field Fp_field={_Fp_red,_Fp_add,_Fp_rmul,_Fp_neg,
    3153             :                                         _Fp_inv,_Fp_equal0,_Fp_s};
    3154             : 
    3155        4557 : const struct bb_field *get_Fp_field(void **E, GEN p)
    3156             : {
    3157        4557 :   *E = (void*)p; return &Fp_field;
    3158             : }
    3159             : 
    3160             : /*********************************************************************/
    3161             : /**                                                                 **/
    3162             : /**               ORDER of INTEGERMOD x  in  (Z/nZ)*                **/
    3163             : /**                                                                 **/
    3164             : /*********************************************************************/
    3165             : ulong
    3166        5124 : Fl_order(ulong a, ulong o, ulong p)
    3167             : {
    3168        5124 :   pari_sp av = avma;
    3169             :   GEN m, P, E;
    3170             :   long i;
    3171        5124 :   if (!o) o = p-1;
    3172        5124 :   m = factoru(o);
    3173        5124 :   P = gel(m,1);
    3174        5124 :   E = gel(m,2);
    3175       12831 :   for (i = lg(P)-1; i; i--)
    3176             :   {
    3177        7707 :     ulong j, l=P[i], e=E[i], t = o / upowuu(l,e), y = Fl_powu(a, t, p);
    3178        7707 :     if (y == 1) o = t;
    3179             :     else {
    3180        5089 :       for (j = 1; j < e; j++) { y = Fl_powu(y, l, p); if (y == 1) break; }
    3181        5089 :       o = t *  upowuu(l, j);
    3182             :     }
    3183             :   }
    3184        5124 :   avma = av; return o;
    3185             : }
    3186             : 
    3187             : /*Find the exact order of a assuming a^o==1*/
    3188             : GEN
    3189       10624 : Fp_order(GEN a, GEN o, GEN p) {
    3190       10624 :   if (lgefint(p) == 3 && typ(o) == t_INT && lgefint(o)==3)
    3191             :   {
    3192          21 :     ulong pp = p[2], oo = o[2];
    3193          21 :     return utoi( Fl_order(umodiu(a, pp), oo, pp) );
    3194             :   }
    3195       10603 :   return gen_order(a, o, (void*)p, &Fp_star);
    3196             : }
    3197             : GEN
    3198          56 : Fp_factored_order(GEN a, GEN o, GEN p)
    3199          56 : { return gen_factored_order(a, o, (void*)p, &Fp_star); }
    3200             : 
    3201             : /* return order of a mod p^e, e > 0, pe = p^e */
    3202             : static GEN
    3203          70 : Zp_order(GEN a, GEN p, long e, GEN pe)
    3204             : {
    3205             :   GEN ap, op;
    3206          70 :   if (equaliu(p, 2))
    3207             :   {
    3208          56 :     if (e == 1) return gen_1;
    3209          56 :     if (e == 2) return mod4(a) == 1? gen_1: gen_2;
    3210          49 :     if (mod4(a) == 1)
    3211          14 :       op = gen_1;
    3212             :     else {
    3213          35 :       op = gen_2;
    3214          35 :       a = Fp_sqr(a, pe);
    3215             :     }
    3216             :   } else {
    3217          14 :     ap = (e == 1)? a: remii(a,p);
    3218          14 :     op = Fp_order(ap, subis(p,1), p);
    3219          14 :     if (e == 1) return op;
    3220           0 :     a = Fp_pow(a, op, pe); /* 1 mod p */
    3221             :   }
    3222          49 :   if (equali1(a)) return op;
    3223           7 :   return mulii(op, powiu(p, e - Z_pval(subis(a,1), p)));
    3224             : }
    3225             : 
    3226             : GEN
    3227          63 : znorder(GEN x, GEN o)
    3228             : {
    3229          63 :   pari_sp av = avma;
    3230             :   GEN b, a;
    3231             : 
    3232          63 :   if (typ(x) != t_INTMOD) pari_err_TYPE("znorder [t_INTMOD expected]",x);
    3233          56 :   b = gel(x,1); a = gel(x,2);
    3234          56 :   if (!equali1(gcdii(a,b))) pari_err_COPRIME("znorder", a,b);
    3235          49 :   if (!o)
    3236             :   {
    3237          35 :     GEN fa = Z_factor(b), P = gel(fa,1), E = gel(fa,2);
    3238          35 :     long i, l = lg(P);
    3239          35 :     o = gen_1;
    3240          70 :     for (i = 1; i < l; i++)
    3241             :     {
    3242          35 :       GEN p = gel(P,i);
    3243          35 :       long e = itos(gel(E,i));
    3244             : 
    3245          35 :       if (l == 2)
    3246          35 :         o = Zp_order(a, p, e, b);
    3247             :       else {
    3248           0 :         GEN pe = powiu(p,e);
    3249           0 :         o = lcmii(o, Zp_order(remii(a,pe), p, e, pe));
    3250             :       }
    3251             :     }
    3252          35 :     return gerepileuptoint(av, o);
    3253             :   }
    3254          14 :   return Fp_order(a, o, b);
    3255             : }
    3256             : GEN
    3257           0 : order(GEN x) { return znorder(x, NULL); }
    3258             : 
    3259             : /*********************************************************************/
    3260             : /**                                                                 **/
    3261             : /**               DISCRETE LOGARITHM  in  (Z/nZ)*                   **/
    3262             : /**                                                                 **/
    3263             : /*********************************************************************/
    3264             : static GEN
    3265       56805 : Fp_log_halfgcd(ulong bnd, GEN C, GEN g, GEN p)
    3266             : {
    3267       56805 :   pari_sp av = avma;
    3268             :   GEN h1, h2, F, G;
    3269       56805 :   if (!Fp_ratlift(g,p,C,shifti(C,-1),&h1,&h2)) return NULL;
    3270       34174 :   if ((F = Z_issmooth_fact(h1, bnd)) && (G = Z_issmooth_fact(h2, bnd)))
    3271             :   {
    3272         231 :     GEN M = cgetg(3, t_MAT);
    3273         231 :     gel(M,1) = vecsmall_concat(gel(F, 1),gel(G, 1));
    3274         231 :     gel(M,2) = vecsmall_concat(gel(F, 2),zv_neg_inplace(gel(G, 2)));
    3275         231 :     return gerepileupto(av, M);
    3276             :   }
    3277       33943 :   avma = av; return NULL;
    3278             : }
    3279             : 
    3280             : static GEN
    3281       56805 : Fp_log_find_rel(GEN b, ulong bnd, GEN C, GEN p, GEN *g, long *e)
    3282             : {
    3283             :   GEN rel;
    3284             :   do
    3285             :   {
    3286       56805 :     (*e)++; *g = Fp_mul(*g, b, p);
    3287       56805 :     rel = Fp_log_halfgcd(bnd, C, *g, p);
    3288       56805 :   } while (!rel);
    3289         231 :   return rel;
    3290             : }
    3291             : 
    3292             : struct Fp_log_rel
    3293             : {
    3294             :   GEN rel;
    3295             :   long *sieve;
    3296             :   ulong prmax;
    3297             :   long nbrel, nbmax;
    3298             : };
    3299             : 
    3300             : /* add u^e */
    3301             : static long
    3302       37401 : addifsmooth1(struct Fp_log_rel *r, GEN h, long u, long e)
    3303             : {
    3304       37401 :   pari_sp av = avma;
    3305             :   GEN z;
    3306       37401 :   if ((z = Z_issmooth_fact(h, r->prmax)))
    3307             :   {
    3308        2485 :     long off = r->prmax+1;
    3309        2485 :     GEN F = cgetg(3, t_MAT);
    3310        2485 :     gel(F,1) = vecsmall_append(gel(z,1), off+u);
    3311        2485 :     gel(F,2) = vecsmall_append(gel(z,2), e);
    3312        2485 :     gel(r->rel,++r->nbrel) = gerepileupto(av, F);
    3313             :   }
    3314       37401 :   return r->nbrel==r->nbmax;
    3315             : }
    3316             : 
    3317             : /* add u^-1 v^-1 */
    3318             : static long
    3319      253358 : addifsmooth2(struct Fp_log_rel *r, GEN h, long u, long v)
    3320             : {
    3321      253358 :   pari_sp av = avma;
    3322             :   GEN z;
    3323      253358 :   if ((z = Z_issmooth_fact(h, r->prmax)))
    3324             :   {
    3325       95011 :     long off = r->prmax+1;
    3326       95011 :     GEN P = mkvecsmall2(off+u,off+v), E = mkvecsmall2(-1,-1);
    3327       95011 :     GEN F = cgetg(3, t_MAT);
    3328       95011 :     gel(F,1) = vecsmall_concat(gel(z,1), P);
    3329       95011 :     gel(F,2) = vecsmall_concat(gel(z,2), E);
    3330       95011 :     gel(r->rel,++r->nbrel) = gerepileupto(av, F);
    3331             :   }
    3332      253358 :   return r->nbrel==r->nbmax;
    3333             : }
    3334             : 
    3335             : /*
    3336             : Let p=C^2+c
    3337             : Solve h = (C+x)*(C+a)-p = 0 [mod l]
    3338             : h= -c+x*(C+a)+C*a = 0  [mod l]
    3339             : x = (c-C*a)/(C+a) [mod l]
    3340             : h = -c+C*(x+a)+a*x
    3341             : */
    3342             : 
    3343             : static void
    3344       36862 : Fp_log_sieve_h(struct Fp_log_rel *r, GEN C, GEN c, GEN Ci, GEN ci, long a, GEN pr, GEN sz)
    3345             : {
    3346       36862 :   long th = expi(C), n = lg(pr)-1;
    3347             :   long i,j;
    3348       36862 :   if (addifsmooth1(r, addis(C,a), a, -1)) return;
    3349    25526508 :   for(j=0; j<=a; j++)
    3350    25489646 :     r->sieve[j]=0;
    3351    16626820 :   for(i=1; i<=n; i++)
    3352             :   {
    3353    16589958 :     ulong li = pr[i], s = sz[i], al = a % li;
    3354    16589958 :     ulong u, iv = Fl_invsafe(Fl_add(Ci[i],al,li),li);
    3355    16589958 :     if (!iv) continue;
    3356    16504691 :     u = Fl_mul(Fl_sub(ci[i],Fl_mul(Ci[i],al,li),li), iv ,li);
    3357    65360344 :     for(j = u; j<=a; j+=li)
    3358    48855653 :       r->sieve[j] += s;
    3359             :   }
    3360       36862 :   th = th - expu(th)-1;
    3361    25465797 :   for(j=0; j<a; j++)
    3362    25428984 :     if (r->sieve[j]>=th)
    3363             :     {
    3364      253358 :       GEN h = addiu(subii(muliu(C,a+j),c), a*j);
    3365      253358 :       if (addifsmooth2(r, h, a, j)) return;
    3366             :     }
    3367             :   /* j = a */
    3368       36813 :     if (r->sieve[a]>=th)
    3369             :     {
    3370         539 :       GEN h = addiu(subii(muliu(C,2*a),c), a*a);
    3371         539 :       if (addifsmooth1(r, h, a, -2)) return;
    3372             :     }
    3373             : }
    3374             : 
    3375             : static GEN
    3376        1367 : _psi(void*E, GEN y)
    3377             : {
    3378        1367 :   GEN lx = (GEN) E;
    3379        1367 :   long prec = realprec(lx);
    3380        1367 :   GEN ly = glog(y, prec);
    3381        1367 :   GEN u = gdiv(lx, ly);
    3382        1367 :   return gsub(gdiv(y ,ly), gpow(u, u, prec));
    3383             : }
    3384             : 
    3385             : static GEN
    3386          49 : opt_param(GEN x, long prec)
    3387             : {
    3388          49 :   return zbrent((void*)glog(x,prec), _psi, gen_2, x, prec);
    3389             : }
    3390             : 
    3391             : static GEN
    3392          49 : check_kernel(long nbg, long N, long prmax, GEN C, GEN M, GEN p, GEN m)
    3393             : {
    3394          49 :   pari_sp av = avma;
    3395          49 :   long lM = lg(M)-1, nbcol = lM;
    3396             :   for (;;)
    3397             :   {
    3398          84 :     GEN K = FpMs_leftkernel_elt_col(M, nbcol, N, m);
    3399          84 :     long i, f=0;
    3400          84 :     long l = lg(K), lm = lgefint(m);
    3401          84 :     GEN idx = diviiexact(subis(p,1),m), g;
    3402             :     pari_timer ti;
    3403          84 :     if (DEBUGLEVEL) timer_start(&ti);
    3404         133 :     for(i=1; i<l; i++)
    3405         133 :       if (signe(gel(K,i)))
    3406          84 :         break;
    3407          84 :     g = Fp_pow(utoi(i), idx, p);
    3408          84 :     K = FpC_Fp_mul(K, Fp_inv(gel(K,i), m), m);
    3409      129185 :     for(i=1; i<l; i++)
    3410             :     {
    3411      129101 :       GEN k = gel(K,i);
    3412      129101 :       GEN j = i<=prmax ? utoi(i): addis(C,i-(prmax+1));
    3413      129101 :       if (signe(k)==0 || !equalii(Fp_pow(g, k, p),
    3414             :             Fp_pow(j, idx, p)))
    3415       80234 :         gel(K,i) = cgetineg(lm);
    3416             :       else
    3417       48867 :         f++;
    3418             :     }
    3419          84 :     if (DEBUGLEVEL) timer_printf(&ti,"found %ld logs", f);
    3420         133 :     if(f > (nbg>>1)) return gerepileupto(av, K);
    3421        9198 :     for(i=1; i<=nbcol; i++)
    3422             :     {
    3423        9163 :       long a = 1+random_Fl(lM);
    3424        9163 :       swap(gel(M,a),gel(M,i));
    3425             :     }
    3426          35 :     if (4*nbcol>5*nbg) nbcol = nbcol*9/10;
    3427          35 :   }
    3428             : }
    3429             : 
    3430             : static GEN
    3431          98 : Fp_log_find_ind(GEN a, GEN K, long prmax, GEN C, GEN p, GEN m)
    3432             : {
    3433          98 :   pari_sp av=avma;
    3434          98 :   GEN aa = gen_1;
    3435          98 :   long AV = 0;
    3436             :   for(;;)
    3437             :   {
    3438         231 :     GEN A = Fp_log_find_rel(a, prmax, C, p, &aa, &AV);
    3439         231 :     GEN F = gel(A,1), E = gel(A,2);
    3440         231 :     GEN Ao = gen_0;
    3441         231 :     long i, l = lg(F);
    3442        1162 :     for(i=1; i<l; i++)
    3443             :     {
    3444        1064 :       GEN Ki = gel(K,F[i]);
    3445        1064 :       if (signe(Ki)<0) break;
    3446         931 :       Ao = addii(Ao, mulis(Ki, E[i]));
    3447             :     }
    3448         329 :     if (i==l) return Fp_div(Ao, utoi(AV), m);
    3449         133 :     aa = gerepileuptoint(av, aa);
    3450         133 :   }
    3451             : }
    3452             : 
    3453             : static GEN
    3454          49 : Fp_log_index(GEN a, GEN b, GEN m, GEN p)
    3455             : {
    3456          49 :   pari_sp av = avma, av2;
    3457             :   long i, nbi, nbrow;
    3458             :   GEN C, c, Ci, ci, pr, sz, l, Ao, Bo, K, d, p_1;
    3459             :   pari_timer ti;
    3460             :   struct Fp_log_rel r;
    3461          49 :   ulong bnds = itou(roundr_safe(opt_param(sqrti(p),DEFAULTPREC)));
    3462          49 :   ulong bnd = 4*bnds;
    3463          49 :   if (!bnds || cmpii(sqru(bnds),m)>=0) return NULL;
    3464             : 
    3465          49 :   p_1 = subiu(p,1);
    3466          49 :   if (!is_pm1(gcdii(m,diviiexact(p_1,m))))
    3467           0 :     m = diviiexact(p_1, coprime_part(p_1, m));
    3468          49 :   pr = primes_upto_zv(bnd);
    3469          49 :   nbi = lg(pr)-1;
    3470          49 :   if (DEBUGLEVEL)
    3471             :   {
    3472           0 :     err_printf("bnd=%lu Size FB=%ld\n", bnd, nbi);
    3473           0 :     timer_start(&ti);
    3474             :   }
    3475          49 :   C = sqrtremi(p, &c);
    3476          49 :   av2 = avma;
    3477          49 :   Ci = cgetg(nbi+1,t_VECSMALL);
    3478          49 :   ci = cgetg(nbi+1,t_VECSMALL);
    3479          49 :   sz = cgetg(nbi+1,t_VECSMALL);
    3480       12236 :   for (i = 1; i <= nbi; ++i)
    3481             :   {
    3482       12187 :     ulong lp = pr[i];
    3483       12187 :     Ci[i] = umodiu(C, lp);
    3484       12187 :     ci[i] = umodiu(c, lp);
    3485       12187 :     sz[i] = expu(lp);
    3486             :   }
    3487          49 :   r.nbrel = 0;
    3488          49 :   r.nbmax = 8*nbi;
    3489          49 :   r.rel = cgetg(r.nbmax+1,t_VEC);
    3490          49 :   r.sieve = cgetg(r.nbmax+2,t_VECSMALL)+1;
    3491          49 :   r.prmax = pr[nbi];
    3492       36911 :   for(i=0; r.nbrel < r.nbmax; i++)
    3493             :   {
    3494       36862 :     Fp_log_sieve_h(&r, C, c, Ci, ci, i, pr, sz);
    3495       36862 :     if (DEBUGLEVEL && (i&127)==0)
    3496           0 :       err_printf("%ld%% ",100*r.nbrel/(r.nbmax));
    3497             :   }
    3498          49 :   nbrow = r.prmax+i;
    3499          49 :   if (DEBUGLEVEL)
    3500             :   {
    3501           0 :     err_printf("\n");
    3502           0 :     timer_printf(&ti," %ld relations, %ld generators", r.nbrel, nbi+i);
    3503             :   }
    3504          49 :   setlg(r.rel,r.nbrel+1);
    3505          49 :   r.rel = gerepileupto(av2, r.rel);
    3506          49 :   K = check_kernel(nbi+nbrow-r.prmax, nbrow, r.prmax, C, r.rel, p, m);
    3507          49 :   if (DEBUGLEVEL) timer_start(&ti);
    3508          49 :   Ao = Fp_log_find_ind(a, K, r.prmax, C, p, m);
    3509          49 :   if (DEBUGLEVEL) timer_printf(&ti," log element");
    3510          49 :   Bo = Fp_log_find_ind(b, K, r.prmax, C, p, m);
    3511          49 :   if (DEBUGLEVEL) timer_printf(&ti," log generator");
    3512          49 :   d = gcdii(Ao,Bo);
    3513          49 :   l = Fp_div(diviiexact(Ao, d) ,diviiexact(Bo, d), m);
    3514          49 :   if (!equalii(a,Fp_pow(b,l,p))) pari_err_BUG("Fp_log_index");
    3515          49 :   return gerepileuptoint(av, l);
    3516             : }
    3517             : 
    3518             : static int
    3519      192262 : Fp_log_use_index(long e, long p)
    3520             : {
    3521      192262 :   return (e >= 27 && 20*(p+6)<=e*e);
    3522             : }
    3523             : 
    3524             : /* Trivial cases a = 1, -1. Return x s.t. g^x = a or [] if no such x exist */
    3525             : static GEN
    3526      218022 : Fp_easylog(void *E, GEN a, GEN g, GEN ord)
    3527             : {
    3528      218022 :   pari_sp av = avma;
    3529      218022 :   GEN p = (GEN)E;
    3530             :   /* assume a reduced mod p, p not necessarily prime */
    3531      218022 :   if (equali1(a)) return gen_0;
    3532             :   /* p > 2 */
    3533      148376 :   if (equalii(subis(p,1), a))  /* -1 */
    3534             :   {
    3535             :     pari_sp av2;
    3536             :     GEN t;
    3537       61716 :     ord = dlog_get_ord(ord);
    3538       61716 :     if (mpodd(ord)) { avma = av; return cgetg(1, t_VEC); } /* no solution */
    3539       61702 :     t = shifti(ord,-1); /* only possible solution */
    3540       61702 :     av2 = avma;
    3541       61702 :     if (!equalii(Fp_pow(g, t, p), a)) { avma = av; return cgetg(1, t_VEC); }
    3542       61674 :     avma = av2; return gerepileuptoint(av, t);
    3543             :   }
    3544       86660 :   if (typ(ord)==t_INT && BPSW_psp(p) && Fp_log_use_index(expi(ord),expi(p)))
    3545          49 :     return Fp_log_index(a, g, ord, p);
    3546       86611 :   avma = av; return NULL; /* not easy */
    3547             : }
    3548             : 
    3549             : GEN
    3550      205479 : Fp_log(GEN a, GEN g, GEN ord, GEN p)
    3551             : {
    3552      205479 :   GEN v = dlog_get_ordfa(ord);
    3553      205451 :   GEN F = gmael(v,2,1);
    3554      205451 :   long lF = lg(F)-1, lmax;
    3555      205451 :   if (lF == 0) return equali1(a)? gen_0: cgetg(1, t_VEC);
    3556      154351 :   lmax = expi(gel(F,lF));
    3557      154351 :   if (BPSW_psp(p) && Fp_log_use_index(lmax,expi(p)))
    3558          49 :     v = mkvec2(gel(v,1),ZM_famat_limit(gel(v,2),int2n(27)));
    3559      154351 :   return gen_PH_log(a,g,v,(void*)p,&Fp_star);
    3560             : }
    3561             : 
    3562             : /* find x such that h = g^x mod N > 1, N = prod_{i <= l} P[i]^E[i], P[i] prime.
    3563             :  * PHI[l] = eulerphi(N / P[l]^E[l]).   Destroys P/E */
    3564             : static GEN
    3565         112 : znlog_rec(GEN h, GEN g, GEN N, GEN P, GEN E, GEN PHI)
    3566             : {
    3567         112 :   long l = lg(P) - 1, e = E[l];
    3568         112 :   GEN p = gel(P, l), phi = gel(PHI,l), pe = e == 1? p: powiu(p, e);
    3569             :   GEN a,b, hp,gp, hpe,gpe, ogpe; /* = order(g mod p^e) | p^(e-1)(p-1) */
    3570             : 
    3571         112 :   if (l == 1) {
    3572          84 :     hpe = h;
    3573          84 :     gpe = g;
    3574             :   } else {
    3575          28 :     hpe = modii(h, pe);
    3576          28 :     gpe = modii(g, pe);
    3577             :   }
    3578         112 :   if (e == 1) {
    3579          28 :     hp = hpe;
    3580          28 :     gp = gpe;
    3581             :   } else {
    3582          84 :     hp = remii(hpe, p);
    3583          84 :     gp = remii(gpe, p);
    3584             :   }
    3585         112 :   if (hp == gen_0 || gp == gen_0) return NULL;
    3586          91 :   if (equaliu(p, 2))
    3587             :   {
    3588          35 :     GEN n = int2n(e);
    3589          35 :     ogpe = Zp_order(gpe, gen_2, e, n);
    3590          35 :     a = Fp_log(hpe, gpe, ogpe, n);
    3591          35 :     if (typ(a) != t_INT) return NULL;
    3592             :   }
    3593             :   else
    3594             :   { /* Avoid black box groups: (Z/p^2)^* / (Z/p)^* ~ (Z/pZ, +), where DL
    3595             :        is trivial */
    3596             :     /* [order(gp), factor(order(gp))] */
    3597          56 :     GEN v = Fp_factored_order(gp, subis(p,1), p);
    3598          56 :     GEN ogp = gel(v,1);
    3599          56 :     if (!equali1(Fp_pow(hp, ogp, p))) return NULL;
    3600          56 :     a = Fp_log(hp, gp, v, p);
    3601          56 :     if (typ(a) != t_INT) return NULL;
    3602          56 :     if (e == 1) ogpe = ogp;
    3603             :     else
    3604             :     { /* find a s.t. g^a = h (mod p^e), p odd prime, e > 0, (h,p) = 1 */
    3605             :       /* use p-adic log: O(log p + e) mul*/
    3606             :       long vpogpe, vpohpe;
    3607             : 
    3608          28 :       hpe = Fp_mul(hpe, Fp_pow(gpe, negi(a), pe), pe);
    3609          28 :       gpe = Fp_pow(gpe, ogp, pe);
    3610             :       /* g,h = 1 mod p; compute b s.t. h = g^b */
    3611             : 
    3612             :       /* v_p(order g mod pe) */
    3613          28 :       vpogpe = equali1(gpe)? 0: e - Z_pval(subis(gpe,1), p);
    3614             :       /* v_p(order h mod pe) */
    3615          28 :       vpohpe = equali1(hpe)? 0: e - Z_pval(subis(hpe,1), p);
    3616          28 :       if (vpohpe > vpogpe) return NULL;
    3617             : 
    3618          28 :       ogpe = mulii(ogp, powiu(p, vpogpe)); /* order g mod p^e */
    3619          28 :       if (is_pm1(gpe)) return is_pm1(hpe)? a: NULL;
    3620          28 :       b = gdiv(Qp_log(cvtop(hpe, p, e)), Qp_log(cvtop(gpe, p, e)));
    3621          28 :       a = addii(a, mulii(ogp, padic_to_Q(b)));
    3622             :     }
    3623             :   }
    3624             :   /* gp^a = hp => x = a mod ogpe => generalized Pohlig-Hellman strategy */
    3625          77 :   if (l == 1) return a;
    3626             : 
    3627          28 :   N = diviiexact(N, pe); /* make N coprime to p */
    3628          28 :   h = Fp_mul(h, Fp_pow(g, modii(negi(a), phi), N), N);
    3629          28 :   g = Fp_pow(g, modii(ogpe, phi), N);
    3630          28 :   setlg(P, l); /* remove last element */
    3631          28 :   setlg(E, l);
    3632          28 :   b = znlog_rec(h, g, N, P, E, PHI);
    3633          28 :   if (!b) return NULL;
    3634          28 :   return addmulii(a, b, ogpe);
    3635             : }
    3636             : 
    3637             : static GEN
    3638          84 : get_PHI(GEN P, GEN E)
    3639             : {
    3640          84 :   long i, l = lg(P);
    3641          84 :   GEN PHI = cgetg(l, t_VEC);
    3642          84 :   gel(PHI,1) = gen_1;
    3643         112 :   for (i=1; i<l-1; i++)
    3644             :   {
    3645          28 :     GEN t, p = gel(P,i);
    3646          28 :     long e = E[i];
    3647          28 :     t = mulii(powiu(p, e-1), subis(p,1));
    3648          28 :     if (i > 1) t = mulii(t, gel(PHI,i));
    3649          28 :     gel(PHI,i+1) = t;
    3650             :   }
    3651          84 :   return PHI;
    3652             : }
    3653             : 
    3654             : GEN
    3655         217 : znlog(GEN h, GEN g, GEN o)
    3656             : {
    3657         217 :   pari_sp av = avma;
    3658             :   GEN N, fa, P, E, x;
    3659         217 :   switch (typ(g))
    3660             :   {
    3661             :     case t_PADIC:
    3662             :     {
    3663          28 :       GEN p = gel(g,2);
    3664          28 :       long v = valp(g);
    3665          28 :       if (v < 0) pari_err_DIM("znlog");
    3666          28 :       if (v > 0) {
    3667           0 :         long k = gvaluation(h, p);
    3668           0 :         if (k % v) return cgetg(1,t_VEC);
    3669           0 :         k /= v;
    3670           0 :         if (!gequal(h, gpowgs(g,k))) { avma = av; return cgetg(1,t_VEC); }
    3671           0 :         avma = av; return stoi(k);
    3672             :       }
    3673          28 :       N = gel(g,3);
    3674          28 :       g = Rg_to_Fp(g, N);
    3675          28 :       break;
    3676             :     }
    3677             :     case t_INTMOD:
    3678         189 :       N = gel(g,1);
    3679         189 :       g = gel(g,2); break;
    3680           0 :     default: pari_err_TYPE("znlog", g);
    3681           0 :       return NULL; /* not reached */
    3682             :   }
    3683         217 :   if (equali1(N)) { avma = av; return gen_0; }
    3684         217 :   h = Rg_to_Fp(h, N);
    3685         210 :   if (o) return gerepileupto(av, Fp_log(h, g, o, N));
    3686          84 :   fa = Z_factor(N);
    3687          84 :   P = gel(fa,1);
    3688          84 :   E = vec_to_vecsmall(gel(fa,2));
    3689          84 :   x = znlog_rec(h, g, N, P, E, get_PHI(P,E));
    3690          84 :   if (!x) { avma = av; return cgetg(1,t_VEC); }
    3691          49 :   return gerepileuptoint(av, x);
    3692             : }
    3693             : 
    3694             : GEN
    3695       62286 : Fp_sqrtn(GEN a, GEN n, GEN p, GEN *zeta)
    3696             : {
    3697       62286 :   a = modii(a,p);
    3698       62286 :   if (!signe(a))
    3699             :   {
    3700       48034 :     if (zeta) *zeta = gen_1;
    3701       48034 :     if (signe(n) < 0) pari_err_INV("Fp_sqrtn", mkintmod(gen_0,p));
    3702       48027 :     return gen_0;
    3703             :   }
    3704       14252 :   if (equaliu(n,2))
    3705             :   {
    3706        1890 :     if (zeta) *zeta = addis(p,-1);
    3707        1890 :     return Fp_sqrt(a,p);
    3708             :   }
    3709       12362 :   return gen_Shanks_sqrtn(a,n,addis(p,-1),zeta,(void*)p,&Fp_star);
    3710             : }
    3711             : 
    3712             : /*********************************************************************/
    3713             : /**                                                                 **/
    3714             : /**                    FUNDAMENTAL DISCRIMINANTS                    **/
    3715             : /**                                                                 **/
    3716             : /*********************************************************************/
    3717             : long
    3718       14371 : isfundamental(GEN x) {
    3719       14371 :   if (typ(x) != t_INT) pari_err_TYPE("isfundamental",x);
    3720       14371 :   return Z_isfundamental(x);
    3721             : }
    3722             : 
    3723             : /* x fundamental ? */
    3724             : long
    3725        7209 : uposisfundamental(ulong x)
    3726             : {
    3727        7209 :   ulong r = x & 15; /* x mod 16 */
    3728        7209 :   if (!r) return 0;
    3729        6817 :   switch(r & 3)
    3730             :   { /* x mod 4 */
    3731        1303 :     case 0: return (r == 4)? 0: uissquarefree(x >> 2);
    3732        2234 :     case 1: return uissquarefree(x);
    3733        3280 :     default: return 0;
    3734             :   }
    3735             : }
    3736             : /* -x fundamental ? */
    3737             : long
    3738       12545 : unegisfundamental(ulong x)
    3739             : {
    3740       12545 :   ulong r = x & 15; /* x mod 16 */
    3741       12545 :   if (!r) return 0;
    3742       12020 :   switch(r & 3)
    3743             :   { /* x mod 4 */
    3744        1982 :     case 0: return (r == 12)? 0: uissquarefree(x >> 2);
    3745        6912 :     case 3: return uissquarefree(x);
    3746        3126 :     default: return 0;
    3747             :   }
    3748             : }
    3749             : long
    3750        1078 : sisfundamental(long x)
    3751        1078 : { return x < 0? unegisfundamental((ulong)(-x)): uposisfundamental(x); }
    3752             : 
    3753             : long
    3754       14931 : Z_isfundamental(GEN x)
    3755             : {
    3756             :   long r;
    3757       14931 :   switch(lgefint(x))
    3758             :   {
    3759           0 :     case 2: return 0;
    3760       32415 :     case 3: return signe(x) < 0? unegisfundamental(x[2])
    3761       19494 :                                : uposisfundamental(x[2]);
    3762             :   }
    3763        2010 :   r = mod16(x);
    3764        2010 :   if (!r) return 0;
    3765        1884 :   if ((r & 3) == 0)
    3766             :   {
    3767             :     pari_sp av;
    3768         376 :     r >>= 2; /* |x|/4 mod 4 */
    3769         376 :     if (signe(x) < 0) r = 4-r;
    3770         376 :     if (r == 1) return 0;
    3771         250 :     av = avma;
    3772         250 :     r = Z_issquarefree( shifti(x,-2) );
    3773         250 :     avma = av; return r;
    3774             :   }
    3775        1508 :   r &= 3; /* |x| mod 4 */
    3776        1508 :   if (signe(x) < 0) r = 4-r;
    3777        1508 :   return (r==1) ? Z_issquarefree(x) : 0;
    3778             : }
    3779             : 
    3780             : GEN
    3781           7 : quaddisc(GEN x)
    3782             : {
    3783           7 :   const pari_sp av = avma;
    3784           7 :   long i,r,tx=typ(x);
    3785             :   GEN P,E,f,s;
    3786             : 
    3787           7 :   if (!is_rational_t(tx)) pari_err_TYPE("quaddisc",x);
    3788           7 :   f = factor(x);
    3789           7 :   P = gel(f,1);
    3790           7 :   E = gel(f,2); s = gen_1;
    3791          35 :   for (i=1; i<lg(P); i++)
    3792          28 :     if (odd(mael(E,i,2))) s = mulii(s,gel(P,i));
    3793           7 :   r = mod4(s); if (gsigne(x) < 0) r = 4-r;
    3794           7 :   if (r>1) s = shifti(s,2);
    3795           7 :   return gerepileuptoint(av, s);
    3796             : }
    3797             : 
    3798             : /*********************************************************************/
    3799             : /**                                                                 **/
    3800             : /**                              FACTORIAL                          **/
    3801             : /**                                                                 **/
    3802             : /*********************************************************************/
    3803             : /* return a * (a+1) * ... * b. Assume a <= b  [ note: factoring out powers of 2
    3804             :  * first is slower ... ] */
    3805             : GEN
    3806      979630 : mulu_interval(ulong a, ulong b)
    3807             : {
    3808      979630 :   pari_sp av = avma;
    3809             :   ulong k, l, N, n;
    3810             :   long lx;
    3811             :   GEN x;
    3812             : 
    3813      979630 :   if (!a) return gen_0;
    3814      979630 :   n = b - a + 1;
    3815      979630 :   if (n < 61)
    3816             :   {
    3817      974603 :     x = utoi(a);
    3818      974606 :     for (k=a+1; k<=b; k++) x = mului(k,x);
    3819      974607 :     return gerepileuptoint(av, x);
    3820             :   }
    3821        5027 :   lx = 1; x = cgetg(2 + n/2, t_VEC);
    3822        5027 :   N = b + a;
    3823      722381 :   for (k = a;; k++)
    3824             :   {
    3825      722381 :     l = N - k; if (l <= k) break;
    3826      717354 :     gel(x,lx++) = muluu(k,l);
    3827      717354 :   }
    3828        5027 :   if (l == k) gel(x,lx++) = utoipos(k);
    3829        5027 :   setlg(x, lx);
    3830        5027 :   return gerepileuptoint(av, ZV_prod(x));
    3831             : }
    3832             : GEN
    3833         245 : muls_interval(long a, long b)
    3834             : {
    3835         245 :   pari_sp av = avma;
    3836         245 :   long lx, k, l, N, n = b - a + 1;
    3837             :   GEN x;
    3838             : 
    3839         245 :   if (a <= 0 && b >= 0) return gen_0;
    3840          98 :   if (n < 61)
    3841             :   {
    3842          98 :     x = stoi(a);
    3843          98 :     for (k=a+1; k<=b; k++) x = mulsi(k,x);
    3844          98 :     return gerepileuptoint(av, x);
    3845             :   }
    3846           0 :   lx = 1; x = cgetg(2 + n/2, t_VEC);
    3847           0 :   N = b + a;
    3848           0 :   for (k = a;; k++)
    3849             :   {
    3850           0 :     l = N - k; if (l <= k) break;
    3851           0 :     gel(x,lx++) = mulss(k,l);
    3852           0 :   }
    3853           0 :   if (l == k) gel(x,lx++) = stoi(k);
    3854           0 :   setlg(x, lx);
    3855           0 :   return gerepileuptoint(av, ZV_prod(x));
    3856             : }
    3857             : 
    3858             : GEN
    3859     2295568 : mpfact(long n)
    3860             : {
    3861     2295568 :   if (n < 2)
    3862             :   {
    3863     1436079 :     if (n < 0) pari_err_DOMAIN("factorial", "argument","<",gen_0,stoi(n));
    3864     1436079 :     return gen_1;
    3865             :   }
    3866      859489 :   return mulu_interval(2UL, (ulong)n);
    3867             : }
    3868             : 
    3869             : /*******************************************************************/
    3870             : /**                                                               **/
    3871             : /**                      LUCAS & FIBONACCI                        **/
    3872             : /**                                                               **/
    3873             : /*******************************************************************/
    3874             : static void
    3875          56 : lucas(ulong n, GEN *a, GEN *b)
    3876             : {
    3877             :   GEN z, t, zt;
    3878         112 :   if (!n) { *a = gen_2; *b = gen_1; return; }
    3879          49 :   lucas(n >> 1, &z, &t); zt = mulii(z, t);
    3880          49 :   switch(n & 3) {
    3881          14 :     case  0: *a = addsi(-2,sqri(z)); *b = addsi(-1,zt); break;
    3882          14 :     case  1: *a = addsi(-1,zt);      *b = addsi(2,sqri(t)); break;
    3883           7 :     case  2: *a = addsi(2,sqri(z));  *b = addsi(1,zt); break;
    3884          14 :     case  3: *a = addsi(1,zt);       *b = addsi(-2,sqri(t));
    3885             :   }
    3886             : }
    3887             : 
    3888             : GEN
    3889           7 : fibo(long n)
    3890             : {
    3891           7 :   pari_sp av = avma;
    3892             :   GEN a, b;
    3893           7 :   if (!n) return gen_0;
    3894           7 :   lucas((ulong)(labs(n)-1), &a, &b);
    3895           7 :   a = diviuexact(addii(shifti(a,1),b), 5);
    3896           7 :   if (n < 0 && !odd(n)) setsigne(a, -1);
    3897           7 :   return gerepileuptoint(av, a);
    3898             : }
    3899             : 
    3900             : /*******************************************************************/
    3901             : /*                                                                 */
    3902             : /*                      CONTINUED FRACTIONS                        */
    3903             : /*                                                                 */
    3904             : /*******************************************************************/
    3905             : static GEN
    3906      359179 : icopy_lg(GEN x, long l)
    3907             : {
    3908      359179 :   long lx = lgefint(x);
    3909             :   GEN y;
    3910             : 
    3911      359179 :   if (lx >= l) return icopy(x);
    3912          35 :   y = cgeti(l); affii(x, y); return y;
    3913             : }
    3914             : 
    3915             : /* continued fraction of a/b. If y != NULL, stop when partial quotients
    3916             :  * differ from y */
    3917             : static GEN
    3918      359468 : Qsfcont(GEN a, GEN b, GEN y, ulong k)
    3919             : {
    3920             :   GEN  z, c;
    3921      359468 :   ulong i, l, ly = lgefint(b);
    3922             : 
    3923             :   /* times 1 / log2( (1+sqrt(5)) / 2 )  */
    3924      359468 :   l = (ulong)(3 + bit_accuracy_mul(ly, 1.44042009041256));
    3925      359468 :   if (k > 0 && k+1 > 0 && l > k+1) l = k+1; /* beware overflow */
    3926      359468 :   if (l > LGBITS) l = LGBITS;
    3927             : 
    3928      359468 :   z = cgetg(l,t_VEC);
    3929      359468 :   l--;
    3930      359468 :   if (y) {
    3931         289 :     pari_sp av = avma;
    3932         289 :     if (l >= (ulong)lg(y)) l = lg(y)-1;
    3933       19467 :     for (i = 1; i <= l; i++)
    3934             :     {
    3935       19289 :       GEN q = gel(y,i);
    3936       19289 :       gel(z,i) = q;
    3937       19289 :       c = b; if (!gequal1(q)) c = mulii(q, b);
    3938       19289 :       c = subii(a, c);
    3939       19289 :       if (signe(c) < 0)
    3940             :       { /* partial quotient too large */
    3941         110 :         c = addii(c, b);
    3942         110 :         if (signe(c) >= 0) i++; /* by 1 */
    3943         110 :         break;
    3944             :       }
    3945       19179 :       if (cmpii(c, b) >= 0)
    3946             :       { /* partial quotient too small */
    3947           1 :         c = subii(c, b);
    3948           1 :         if (cmpii(c, b) < 0) {
    3949             :           /* by 1. If next quotient is 1 in y, add 1 */
    3950           0 :           if (i < l && equali1(gel(y,i+1))) gel(z,i) = addis(q,1);
    3951           0 :           i++;
    3952             :         }
    3953           1 :         break;
    3954             :       }
    3955       19178 :       if ((i & 0xff) == 0) gerepileall(av, 2, &b, &c);
    3956       19178 :       a = b; b = c;
    3957             :     }
    3958             :   } else {
    3959      359179 :     a = icopy_lg(a, ly);
    3960      359179 :     b = icopy(b);
    3961     1254340 :     for (i = 1; i <= l; i++)
    3962             :     {
    3963     1254070 :       gel(z,i) = truedvmdii(a,b,&c);
    3964     1254070 :       if (c == gen_0) { i++; break; }
    3965      895161 :       affii(c, a); cgiv(c); c = a;
    3966      895161 :       a = b; b = c;
    3967             :     }
    3968             :   }
    3969      359468 :   i--;
    3970      359468 :   if (i > 1 && gequal1(gel(z,i)))
    3971             :   {
    3972          79 :     cgiv(gel(z,i)); --i;
    3973          79 :     gel(z,i) = addsi(1, gel(z,i)); /* unclean: leave old z[i] on stack */
    3974             :   }
    3975      359468 :   setlg(z,i+1); return z;
    3976             : }
    3977             : 
    3978             : static GEN
    3979           0 : sersfcont(GEN a, GEN b, long k)
    3980             : {
    3981           0 :   long i, l = typ(a) == t_POL? lg(a): 3;
    3982             :   GEN y, c;
    3983           0 :   if (lg(b) > l) l = lg(b);
    3984           0 :   if (k > 0 && l > k+1) l = k+1;
    3985           0 :   y = cgetg(l,t_VEC);
    3986           0 :   for (i=1; i<l; i++)
    3987             :   {
    3988           0 :     gel(y,i) = poldivrem(a,b,&c);
    3989           0 :     if (gequal0(c)) { i++; break; }
    3990           0 :     a = b; b = c;
    3991             :   }
    3992           0 :   setlg(y, i); return y;
    3993             : }
    3994             : 
    3995             : GEN
    3996      359473 : gboundcf(GEN x, long k)
    3997             : {
    3998             :   pari_sp av;
    3999      359473 :   long tx = typ(x), e;
    4000             :   GEN y, a, b, c;
    4001             : 
    4002      359473 :   if (k < 0) pari_err_DOMAIN("gboundcf","nmax","<",gen_0,stoi(k));
    4003      359466 :   if (is_scalar_t(tx))
    4004             :   {
    4005      359466 :     if (gequal0(x)) return mkvec(gen_0);
    4006      359186 :     switch(tx)
    4007             :     {
    4008           0 :       case t_INT: return mkveccopy(x);
    4009             :       case t_REAL:
    4010         296 :         av = avma;
    4011         296 :         c = mantissa_real(x,&e);
    4012         296 :         if (e < 0) pari_err_PREC("gboundcf");
    4013         289 :         y = int2n(e);
    4014         289 :         a = Qsfcont(c,y, NULL, k);
    4015         289 :         b = addsi(signe(x), c);
    4016         289 :         return gerepilecopy(av, Qsfcont(b,y, a, k));
    4017             : 
    4018             :       case t_FRAC:
    4019      358890 :         av = avma;
    4020      358890 :         return gerepileupto(av, Qsfcont(gel(x,1),gel(x,2), NULL, k));
    4021             :     }
    4022           0 :     pari_err_TYPE("gboundcf",x);
    4023             :   }
    4024             : 
    4025           0 :   switch(tx)
    4026             :   {
    4027           0 :     case t_POL: return mkveccopy(x);
    4028             :     case t_SER:
    4029           0 :       av = avma;
    4030           0 :       return gerepileupto(av, gboundcf(ser2rfrac_i(x), k));
    4031             :     case t_RFRAC:
    4032           0 :       av = avma;
    4033           0 :       return gerepilecopy(av, sersfcont(gel(x,1), gel(x,2), k));
    4034             :   }
    4035           0 :   pari_err_TYPE("gboundcf",x);
    4036           0 :   return NULL; /* not reached */
    4037             : }
    4038             : 
    4039             : static GEN
    4040          14 : sfcont2(GEN b, GEN x, long k)
    4041             : {
    4042          14 :   pari_sp av = avma;
    4043          14 :   long lb = lg(b), tx = typ(x), i;
    4044             :   GEN y,p1;
    4045             : 
    4046          14 :   if (k)
    4047             :   {
    4048           7 :     if (k >= lb) pari_err_DIM("contfrac [too few denominators]");
    4049           0 :     lb = k+1;
    4050             :   }
    4051           7 :   y = cgetg(lb,t_VEC);
    4052           7 :   if (lb==1) return y;
    4053           7 :   if (is_scalar_t(tx))
    4054             :   {
    4055           7 :     if (!is_intreal_t(tx) && tx != t_FRAC) pari_err_TYPE("sfcont2",x);
    4056             :   }
    4057           0 :   else if (tx == t_SER) x = ser2rfrac_i(x);
    4058             : 
    4059           7 :   if (!gequal1(gel(b,1))) x = gmul(gel(b,1),x);
    4060           7 :   for (i = 1;;)
    4061             :   {
    4062          35 :     if (tx == t_REAL)
    4063             :     {
    4064          35 :       long e = expo(x);
    4065          35 :       if (e > 0 && nbits2prec(e+1) > realprec(x)) break;
    4066          35 :       gel(y,i) = floorr(x);
    4067          35 :       p1 = subri(x, gel(y,i));
    4068             :     }
    4069             :     else
    4070             :     {
    4071           0 :       gel(y,i) = gfloor(x);
    4072           0 :       p1 = gsub(x, gel(y,i));
    4073             :     }
    4074          35 :     if (++i >= lb) break;
    4075          28 :     if (gequal0(p1)) break;
    4076          28 :     x = gdiv(gel(b,i),p1);
    4077          28 :   }
    4078           7 :   setlg(y,i);
    4079           7 :   return gerepilecopy(av,y);
    4080             : }
    4081             : 
    4082             : 
    4083             : GEN
    4084           0 : gcf(GEN x) { return gboundcf(x,0); }
    4085             : GEN
    4086           0 : gcf2(GEN b, GEN x) { return contfrac0(x,b,0); }
    4087             : GEN
    4088          49 : contfrac0(GEN x, GEN b, long nmax)
    4089             : {
    4090             :   long tb;
    4091             : 
    4092          49 :   if (!b) return gboundcf(x,nmax);
    4093          28 :   tb = typ(b);
    4094          28 :   if (tb == t_INT) return gboundcf(x,itos(b));
    4095          21 :   if (! is_vec_t(tb)) pari_err_TYPE("contfrac0",b);
    4096          21 :   if (nmax < 0) pari_err_DOMAIN("contfrac","nmax","<",gen_0,stoi(nmax));
    4097          14 :   return sfcont2(b,x,nmax);
    4098             : }
    4099             : 
    4100             : GEN
    4101         126 : contfracpnqn(GEN x, long n)
    4102             : {
    4103         126 :   pari_sp av = avma;
    4104         126 :   long i, lx = lg(x);
    4105             :   GEN M,A,B, p0,p1, q0,q1;
    4106             : 
    4107         126 :   if (lx == 1)
    4108             :   {
    4109          28 :     if (! is_matvec_t(typ(x))) pari_err_TYPE("pnqn",x);
    4110          21 :     if (n >= 0) return cgetg(1,t_MAT);
    4111           7 :     return matid(2);
    4112             :   }
    4113          98 :   switch(typ(x))
    4114             :   {
    4115          56 :     case t_VEC: case t_COL: A = x; B = NULL; break;
    4116             :     case t_MAT:
    4117          42 :       switch(lgcols(x))
    4118             :       {
    4119           0 :         case 2: A = row(x,1); B = NULL; break;
    4120          35 :         case 3: A = row(x,2); B = row(x,1); break;
    4121           7 :         default: pari_err_DIM("pnqn [ nbrows != 1,2 ]");
    4122           0 :                  return NULL; /*not reached*/
    4123             :       }
    4124          35 :       break;
    4125           0 :     default: pari_err_TYPE("pnqn",x);
    4126           0 :       return NULL; /*not reached*/
    4127             :   }
    4128          91 :   p1 = gel(A,1);
    4129          91 :   q1 = B? gel(B,1): gen_1; /* p[0], q[0] */
    4130          91 :   if (n >= 0)
    4131             :   {
    4132          56 :     lx = minss(lx, n+2);
    4133          56 :     if (lx == 2) return gerepilecopy(av, mkmat(mkcol2(p1,q1)));
    4134             :   }
    4135          35 :   else if (lx == 2)
    4136           7 :     return gerepilecopy(av, mkmat2(mkcol2(p1,q1), mkcol2(gen_1,gen_0)));
    4137             :   /* lx >= 3 */
    4138          56 :   p0 = gen_1;
    4139          56 :   q0 = gen_0; /* p[-1], q[-1] */
    4140          56 :   M = cgetg(lx, t_MAT);
    4141          56 :   gel(M,1) = mkcol2(p1,q1);
    4142         217 :   for (i=2; i<lx; i++)
    4143             :   {
    4144         161 :     GEN a = gel(A,i), p2,q2;
    4145         161 :     if (B) {
    4146          84 :       GEN b = gel(B,i);
    4147          84 :       p0 = gmul(b,p0);
    4148          84 :       q0 = gmul(b,q0);
    4149             :     }
    4150         161 :     p2 = gadd(gmul(a,p1),p0); p0=p1; p1=p2;
    4151         161 :     q2 = gadd(gmul(a,q1),q0); q0=q1; q1=q2;
    4152         161 :     gel(M,i) = mkcol2(p1,q1);
    4153             :   }
    4154          56 :   if (n < 0) M = mkmat2(gel(M,lx-1), gel(M,lx-2));
    4155          56 :   return gerepilecopy(av, M);
    4156             : }
    4157             : GEN
    4158           0 : pnqn(GEN x) { return contfracpnqn(x,-1); }
    4159             : /* x = [a0, ..., an] from gboundcf, n >= 0;
    4160             :  * return [[p0, ..., pn], [q0,...,qn]] */
    4161             : GEN
    4162      352730 : ZV_allpnqn(GEN x)
    4163             : {
    4164      352730 :   long i, lx = lg(x);
    4165      352730 :   GEN p0, p1, q0, q1, p2, q2, P,Q, v = cgetg(3,t_VEC);
    4166             : 
    4167      352730 :   gel(v,1) = P = cgetg(lx, t_VEC);
    4168      352730 :   gel(v,2) = Q = cgetg(lx, t_VEC);
    4169      352730 :   p0 = gen_1; q0 = gen_0;
    4170      352730 :   gel(P, 1) = p1 = gel(x,1); gel(Q, 1) = q1 = gen_1;
    4171     1204644 :   for (i=2; i<lx; i++)
    4172             :   {
    4173      851914 :     GEN a = gel(x,i);
    4174      851914 :     gel(P, i) = p2 = addmulii(p0, a, p1); p0 = p1; p1 = p2;
    4175      851914 :     gel(Q, i) = q2 = addmulii(q0, a, q1); q0 = q1; q1 = q2;
    4176             :   }
    4177      352730 :   return v;
    4178             : }
    4179             : 
    4180             : /* write Mod(x,N) as a/b, gcd(a,b) = 1, b <= B (no condition if B = NULL) */
    4181             : static GEN
    4182          35 : mod_to_frac(GEN x, GEN N, GEN B)
    4183             : {
    4184             :   GEN a, b, A;
    4185          35 :   if (B)
    4186             :   {
    4187          21 :     A = divii(shifti(N, -1), B);
    4188             :     /* denominator bound useless, don't use it */
    4189          21 :     if (cmpii(A, B) < 0) B = NULL;
    4190             :   }
    4191          35 :   if (!B)
    4192             :   {
    4193          14 :     A = sqrti(shifti(N, -1));
    4194          14 :     B = A;
    4195             :   }
    4196          35 :   if (!Fp_ratlift(x, N, A,B,&a,&b) || !equali1( gcdii(a,b) )) return NULL;
    4197          21 :   return equali1(b)? a: mkfrac(a,b);
    4198             : }
    4199             : 
    4200             : static GEN
    4201          56 : mod_to_rfrac(GEN x, GEN N, long B)
    4202             : {
    4203             :   GEN a, b;
    4204          56 :   long A, d = degpol(N);
    4205          56 :   if (B >= 0)
    4206             :   {
    4207          21 :     A = d-1 - B;
    4208             :     /* denominator bound useless, don't use it */
    4209          21 :     if (A < B) B = -1;
    4210             :   }
    4211          56 :   if (B < 0)
    4212             :   {
    4213          42 :     B = d >> 1;
    4214          42 :     A = odd(d)? B : B-1;
    4215             :   }
    4216          56 :   if (varn(N) != varn(x)) x = scalarpol(x, varn(N));
    4217          56 :   if (! RgXQ_ratlift(x, N, A, B, &a,&b)) return NULL;
    4218          56 :   if (degpol(RgX_gcd(a,b)) > 0) return NULL;
    4219          49 :   return gdiv(a,b);
    4220             : }
    4221             : 
    4222             : /* k > 0 t_INT, x a t_FRAC, returns the convergent a/b
    4223             :  * of the continued fraction of x with b <= k maximal */
    4224             : static GEN
    4225           0 : bestappr_frac(GEN x, GEN k)
    4226             : {
    4227             :   pari_sp av;
    4228             :   GEN p0, p1, p, q0, q1, q, a, y;
    4229             : 
    4230           0 :   if (cmpii(gel(x,2),k) <= 0) return gcopy(x);
    4231           0 :   av = avma; y = x;
    4232           0 :   p1 = gen_1; p0 = truedvmdii(gel(x,1), gel(x,2), &a); /* = floor(x) */
    4233           0 :   q1 = gen_0; q0 = gen_1;
    4234           0 :   x = mkfrac(a, gel(x,2)); /* = frac(x); now 0<= x < 1 */
    4235             :   for(;;)
    4236             :   {
    4237           0 :     x = ginv(x); /* > 1 */
    4238           0 :     a = typ(x)==t_INT? x: divii(gel(x,1), gel(x,2));
    4239           0 :     if (cmpii(a,k) > 0)
    4240             :     { /* next partial quotient will overflow limits */
    4241             :       GEN n, d;
    4242           0 :       a = divii(subii(k, q1), q0);
    4243           0 :       p = addii(mulii(a,p0), p1); p1=p0; p0=p;
    4244           0 :       q = addii(mulii(a,q0), q1); q1=q0; q0=q;
    4245             :       /* compare |y-p0/q0|, |y-p1/q1| */
    4246           0 :       n = gel(y,1);
    4247           0 :       d = gel(y,2);
    4248           0 :       if (absi_cmp(mulii(q1, subii(mulii(q0,n), mulii(d,p0))),
    4249             :                    mulii(q0, subii(mulii(q1,n), mulii(d,p1)))) < 0)
    4250           0 :                    { p1 = p0; q1 = q0; }
    4251           0 :       break;
    4252             :     }
    4253           0 :     p = addii(mulii(a,p0), p1); p1=p0; p0=p;
    4254           0 :     q = addii(mulii(a,q0), q1); q1=q0; q0=q;
    4255             : 
    4256           0 :     if (cmpii(q0,k) > 0) break;
    4257           0 :     x = gsub(x,a); /* 0 <= x < 1 */
    4258           0 :     if (typ(x) == t_INT) { p1 = p0; q1 = q0; break; } /* x = 0 */
    4259             : 
    4260           0 :   }
    4261           0 :   return gerepileupto(av, gdiv(p1,q1));
    4262             : }
    4263             : /* bestappr(t_REAL != 0), to maximal accuracy */
    4264             : static GEN
    4265          21 : bestappr_real_max(GEN x)
    4266             : {
    4267          21 :   pari_sp av = avma;
    4268             :   GEN p0, p1, p, q0, q1, q, a;
    4269             :   long e;
    4270          21 :   p1 = gen_1; a = p0 = floorr(x); q1 = gen_0; q0 = gen_1;
    4271          21 :   x = subri(x,a); /* 0 <= x < 1 */
    4272          21 :   e = bit_prec(x) - expo(x);
    4273             :   for(;;)
    4274             :   {
    4275             :     long d;
    4276          42 :     if (!signe(x) || expi(q0) > e) { p1 = p0; q1 = q0; break; }
    4277          22 :     x = invr(x); /* > 1 */
    4278          22 :     d = nbits2prec(expo(x) + 1);
    4279          22 :     if (d > lg(x)) { p1 = p0; q1 = q0; break; } /* original x was ~ 0 */
    4280             : 
    4281          21 :     a = truncr(x); /* truncr(x) will NOT raise e_PREC */
    4282          21 :     p = addii(mulii(a,p0), p1); p1=p0; p0=p;
    4283          21 :     q = addii(mulii(a,q0), q1); q1=q0; q0=q;
    4284          21 :     x = subri(x,a); /* 0 <= x < 1 */
    4285          21 :   }
    4286          21 :   return gerepileupto(av, gdiv(p1,q1));
    4287             : }
    4288             : /* k > 0 t_INT, x != 0 a t_REAL, returns the convergent a/b
    4289             :  * of the continued fraction of x with b <= k maximal */
    4290             : static GEN
    4291      195568 : bestappr_real(GEN x, GEN k)
    4292             : {
    4293      195568 :   pari_sp av = avma;
    4294             :   GEN kr, p0, p1, p, q0, q1, q, a, y;
    4295             : 
    4296      195568 :   y = x;
    4297      195568 :   p1 = gen_1; a = p0 = floorr(x);
    4298      195568 :   q1 = gen_0; q0 = gen_1;
    4299      195568 :   x = subri(x,a); /* 0 <= x < 1 */
    4300      195568 :   if (!signe(x)) { cgiv(x); return a; }
    4301      191651 :   kr = itor(k, realprec(x));
    4302             :   for(;;)
    4303             :   {
    4304             :     long d;
    4305      544418 :     x = invr(x); /* > 1 */
    4306      544418 :     if (cmprr(x,kr) > 0)
    4307             :     { /* next partial quotient will overflow limits */
    4308      187197 :       a = divii(subii(k, q1), q0);
    4309      187197 :       p = addii(mulii(a,p0), p1); p1=p0; p0=p;
    4310      187197 :       q = addii(mulii(a,q0), q1); q1=q0; q0=q;
    4311             :       /* compare |y-p0/q0|, |y-p1/q1| */
    4312      187197 :       if (absr_cmp(mulir(q1, subri(mulir(q0,y), p0)),
    4313             :                    mulir(q0, subri(mulir(q1,y), p1))) < 0)
    4314        5059 :                    { p1 = p0; q1 = q0; }
    4315      187197 :       break;
    4316             :     }
    4317      357221 :     d = nbits2prec(expo(x) + 1);
    4318      357221 :     if (d > lg(x)) { p1 = p0; q1 = q0; break; } /* original x was ~ 0 */
    4319             : 
    4320      357218 :     a = truncr(x); /* truncr(x) will NOT raise e_PREC */
    4321      357218 :     p = addii(mulii(a,p0), p1); p1=p0; p0=p;
    4322      357218 :     q = addii(mulii(a,q0), q1); q1=q0; q0=q;
    4323             : 
    4324      357218 :     if (cmpii(q0,k) > 0) break;
    4325      353134 :     x = subri(x,a); /* 0 <= x < 1 */
    4326      353134 :     if (!signe(x)) { p1 = p0; q1 = q0; break; }
    4327      352767 :   }
    4328      191651 :   return gerepileupto(av, gdiv(p1,q1));
    4329             : }
    4330             : 
    4331             : /* k t_INT or NULL */
    4332             : static GEN
    4333      290592 : bestappr_Q(GEN x, GEN k)
    4334             : {
    4335      290592 :   long lx, tx = typ(x), i;
    4336             :   GEN a, y;
    4337             : 
    4338      290592 :   switch(tx)
    4339             :   {
    4340           0 :     case t_INT: return icopy(x);
    4341           0 :     case t_FRAC: return k? bestappr_frac(x, k): gcopy(x);
    4342             :     case t_REAL:
    4343      223051 :       if (!signe(x)) return gen_0;
    4344      195589 :       return k? bestappr_real(x, k): bestappr_real_max(x);
    4345             : 
    4346             :     case t_INTMOD: {
    4347          21 :       pari_sp av = avma;
    4348          21 :       a = mod_to_frac(gel(x,2), gel(x,1), k); if (!a) return NULL;
    4349          14 :       return gerepilecopy(av, a);
    4350             :     }
    4351             :     case t_PADIC: {
    4352          14 :       pari_sp av = avma;
    4353          14 :       long v = valp(x);
    4354          14 :       a = mod_to_frac(gel(x,4), gel(x,3), k); if (!a) return NULL;
    4355           7 :       if (v) a = gmul(a, powis(gel(x,2), v));
    4356           7 :       return gerepilecopy(av, a);
    4357             :     }
    4358             : 
    4359             :     case t_SER:
    4360           0 :       if (ser_isexactzero(x)) return gcopy(x);
    4361             :       /* fall through */
    4362             :     case t_COMPLEX: case t_POLMOD: case t_POL: case t_RFRAC:
    4363             :     case t_VEC: case t_COL: case t_MAT:
    4364       67506 :       y = cgetg_copy(x, &lx);
    4365       67506 :       if (lontyp[tx] == 1) i = 1; else { y[1] = x[1]; i = 2; }
    4366      350109 :       for (; i<lx; i++)
    4367             :       {
    4368      282603 :         a = bestappr_Q(gel(x,i),k); if (!a) return NULL;
    4369      282603 :         gel(y,i) = a;
    4370             :       }
    4371       67506 :       if (tx == t_POL) return normalizepol(y);
    4372       67506 :       if (tx == t_SER) return normalize(y);
    4373       67506 :       return y;
    4374             :   }
    4375           0 :   pari_err_TYPE("bestappr_Q",x);
    4376           0 :   return NULL; /* not reached */
    4377             : }
    4378             : 
    4379             : static GEN
    4380          49 : bestappr_ser(GEN x, long B)
    4381             : {
    4382          49 :   long v = valp(x), lx = lg(x);
    4383             :   GEN N, t;
    4384          49 :   x = normalizepol(ser2pol_i(x, lx));
    4385          49 :   N = monomial(gen_1, lx-2, varn(x));
    4386          49 :   t = mod_to_rfrac(x, N, B); if (!t) return NULL;
    4387          42 :   if (v)
    4388             :   {
    4389             :     GEN a, b;
    4390             :     long vx;
    4391          14 :     if (typ(t) == t_POL) return RgX_mulXn(t, v);
    4392             :     /* t_RFRAC */
    4393          14 :     vx = varn(x);
    4394          14 :     a = gel(t,1);
    4395          14 :     b = gel(t,2);
    4396          14 :     v -= RgX_valrem(b, &b);
    4397          14 :     if (typ(a) == t_POL && varn(a) == vx) v += RgX_valrem(a, &a);
    4398          14 :     if (v < 0) b = RgX_shift(b, -v);
    4399           7 :     else if (v > 0) {
    4400           7 :       if (typ(a) != t_POL || varn(a) != vx) a = scalarpol_shallow(a, vx);
    4401           7 :       a = RgX_shift(a, v);
    4402             :     }
    4403          14 :     t = mkrfraccopy(a, b);
    4404             :   }
    4405          42 :   return t;
    4406             : }
    4407             : static GEN bestappr_RgX(GEN x, long B);
    4408             : /* x t_POLMOD, B >= 0 or < 0 [omit condition on B].
    4409             :  * Look for coprime t_POL a,b, deg(b)<=B, such that a/b = x */
    4410             : static GEN
    4411          63 : bestappr_RgX(GEN x, long B)
    4412             : {
    4413          63 :   long i, lx, tx = typ(x);
    4414             :   GEN y, t;
    4415          63 :   switch(tx)
    4416             :   {
    4417             :     case t_INT: case t_REAL: case t_INTMOD: case t_FRAC:
    4418             :     case t_COMPLEX: case t_PADIC: case t_QUAD: case t_POL:
    4419           0 :       return gcopy(x);
    4420             : 
    4421             :     case t_RFRAC: {
    4422          14 :       pari_sp av = avma;
    4423          14 :       if (B < 0 || degpol(gel(x,2)) <= B) return gcopy(x);
    4424           7 :       x = rfractoser(x, varn(gel(x,2)), 2*B+1);
    4425           7 :       t = bestappr_ser(x, B); if (!t) return NULL;
    4426           7 :       return gerepileupto(av, t);
    4427             :     }
    4428             :     case t_POLMOD: {
    4429           7 :       pari_sp av = avma;
    4430           7 :       t = mod_to_rfrac(gel(x,2), gel(x,1), B); if (!t) return NULL;
    4431           7 :       return gerepileupto(av, t);
    4432             :     }
    4433             :     case t_SER: {
    4434          42 :       pari_sp av = avma;
    4435          42 :       t = bestappr_ser(x, B); if (!t) return NULL;
    4436          35 :       return gerepileupto(av, t);
    4437             :     }
    4438             : 
    4439             :     case t_VEC: case t_COL: case t_MAT:
    4440           0 :       y = cgetg_copy(x, &lx);
    4441           0 :       if (lontyp[tx] == 1) i = 1; else { y[1] = x[1]; i = 2; }
    4442           0 :       for (; i<lx; i++)
    4443             :       {
    4444           0 :         t = bestappr_RgX(gel(x,i),B); if (!t) return NULL;
    4445           0 :         gel(y,i) = t;
    4446             :       }
    4447           0 :       return y;
    4448             :   }
    4449           0 :   pari_err_TYPE("bestappr_RgX",x);
    4450           0 :   return NULL; /* not reached */
    4451             : }
    4452             : 
    4453             : /* allow k = NULL: maximal accuracy */
    4454             : GEN
    4455        7989 : bestappr(GEN x, GEN k)
    4456             : {
    4457        7989 :   pari_sp av = avma;
    4458        7989 :   if (k) { /* replace by floor(k) */
    4459        7954 :     switch(typ(k))
    4460             :     {
    4461             :       case t_INT:
    4462         161 :         break;
    4463             :       case t_REAL: case t_FRAC:
    4464        7793 :         k = floor_safe(k); /* left on stack for efficiency */
    4465        7793 :         if (!signe(k)) k = gen_1;
    4466        7793 :         break;
    4467             :       default:
    4468           0 :         pari_err_TYPE("bestappr [bound type]", k);
    4469           0 :         break;
    4470             :     }
    4471             :   }
    4472        7989 :   x = bestappr_Q(x, k);
    4473        7989 :   if (!x) { avma = av; return cgetg(1,t_VEC); }
    4474        7975 :   return x;
    4475             : }
    4476             : GEN
    4477          63 : bestapprPade(GEN x, long B)
    4478             : {
    4479          63 :   pari_sp av = avma;
    4480          63 :   GEN t = bestappr_RgX(x, B);
    4481          63 :   if (!t) { avma = av; return cgetg(1,t_VEC); }
    4482          56 :   return t;
    4483             : }
    4484             : 
    4485             : /***********************************************************************/
    4486             : /**                                                                   **/
    4487             : /**         FUNDAMENTAL UNIT AND REGULATOR (QUADRATIC FIELDS)         **/
    4488             : /**                                                                   **/
    4489             : /***********************************************************************/
    4490             : 
    4491             : static GEN
    4492          14 : get_quad(GEN f, GEN pol, long r)
    4493             : {
    4494          14 :   GEN p1 = gcoeff(f,1,2), q1 = gcoeff(f,2,2);
    4495          14 :   return mkquad(pol, r? subii(p1,q1): p1, q1);
    4496             : }
    4497             : 
    4498             : /* replace f by f * [a,1; 1,0] */
    4499             : static void
    4500          14 : update_f(GEN f, GEN a)
    4501             : {
    4502             :   GEN p1;
    4503          14 :   p1 = gcoeff(f,1,1);
    4504          14 :   gcoeff(f,1,1) = addii(mulii(a,p1), gcoeff(f,1,2));
    4505          14 :   gcoeff(f,1,2) = p1;
    4506             : 
    4507          14 :   p1 = gcoeff(f,2,1);
    4508          14 :   gcoeff(f,2,1) = addii(mulii(a,p1), gcoeff(f,2,2));
    4509          14 :   gcoeff(f,2,2) = p1;
    4510          14 : }
    4511             : 
    4512             : GEN
    4513           7 : quadunit(GEN x)
    4514             : {
    4515           7 :   pari_sp av = avma, av2;
    4516             :   GEN pol, y, a, u, v, sqd, f;
    4517             :   long r;
    4518             : 
    4519           7 :   check_quaddisc_real(x, &r, "quadunit");
    4520           7 :   pol = quadpoly(x);
    4521           7 :   sqd = sqrti(x); av2 = avma;
    4522           7 :   a = shifti(addsi(r,sqd),-1);
    4523           7 :   f = mkmat2(mkcol2(a, gen_1), mkcol2(gen_1, gen_0)); /* [a,0; 1,0] */
    4524           7 :   u = stoi(r); v = gen_2;
    4525             :   for(;;)
    4526             :   {
    4527             :     GEN u1, v1;
    4528          14 :     u1 = subii(mulii(a,v),u);
    4529          14 :     v1 = divii(subii(x,sqri(u1)),v);
    4530          14 :     if ( equalii(v,v1) ) {
    4531           7 :       y = get_quad(f,pol,r);
    4532           7 :       update_f(f,a);
    4533           7 :       y = gdiv(get_quad(f,pol,r), gconj(y));
    4534           7 :       break;
    4535             :     }
    4536           7 :     a = divii(addii(sqd,u1), v1);
    4537           7 :     if ( equalii(u,u1) ) {
    4538           0 :       y = get_quad(f,pol,r);
    4539           0 :       y = gdiv(y, gconj(y));
    4540           0 :       break;
    4541             :     }
    4542           7 :     update_f(f,a);
    4543           7 :     u = u1; v = v1;
    4544           7 :     if (gc_needed(av2,2))
    4545             :     {
    4546           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"quadunit");
    4547           0 :       gerepileall(av2,4, &a,&f,&u,&v);
    4548             :     }
    4549           7 :   }
    4550           7 :   if (signe(gel(y,3)) < 0) y = gneg(y);
    4551           7 :   return gerepileupto(av, y);
    4552             : }
    4553             : 
    4554             : GEN
    4555          21 : quadregulator(GEN x, long prec)
    4556             : {
    4557          21 :   pari_sp av = avma, av2;
    4558             :   GEN R, rsqd, u, v, sqd;
    4559             :   long r, Rexpo;
    4560             : 
    4561          21 :   check_quaddisc_real(x, &r, "quadregulator");
    4562          21 :   sqd = sqrti(x);
    4563          21 :   rsqd = gsqrt(x,prec);
    4564          21 :   Rexpo = 0; R = real2n(1, prec); /* = 2 */
    4565          21 :   av2 = avma;
    4566          21 :   u = stoi(r); v = gen_2;
    4567             :   for(;;)
    4568             :   {
    4569          70 :     GEN u1 = subii(mulii(divii(addii(u,sqd),v), v), u);
    4570          70 :     GEN v1 = divii(subii(x,sqri(u1)),v);
    4571          70 :     if (equalii(v,v1))
    4572             :     {
    4573           7 :       R = sqrr(R); shiftr_inplace(R, -1);
    4574           7 :       R = mulrr(R, divri(addir(u1,rsqd),v));
    4575           7 :       break;
    4576             :     }
    4577          63 :     if (equalii(u,u1))
    4578             :     {
    4579          14 :       R = sqrr(R); shiftr_inplace(R, -1);
    4580          14 :       break;
    4581             :     }
    4582          49 :     R = mulrr(R, divri(addir(u1,rsqd),v));
    4583          49 :     Rexpo += expo(R); setexpo(R,0);
    4584          49 :     u = u1; v = v1;
    4585          49 :     if (Rexpo & ~EXPOBITS) pari_err_OVERFLOW("quadregulator [exponent]");
    4586          49 :     if (gc_needed(av2,2))
    4587             :     {
    4588           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"quadregulator");
    4589           0 :       gerepileall(av2,3, &R,&u,&v);
    4590             :     }
    4591          49 :   }
    4592          21 :   R = logr_abs(divri(R,v));
    4593          21 :   if (Rexpo)
    4594             :   {
    4595          21 :     GEN t = mulsr(Rexpo, mplog2(prec));
    4596          21 :     shiftr_inplace(t, 1);
    4597          21 :     R = addrr(R,t);
    4598             :   }
    4599          21 :   return gerepileuptoleaf(av, R);
    4600             : }
    4601             : 
    4602             : /*************************************************************************/
    4603             : /**                                                                     **/
    4604             : /**                            CLASS NUMBER                             **/
    4605             : /**                                                                     **/
    4606             : /*************************************************************************/
    4607             : 
    4608             : int
    4609    10628904 : qfb_equal1(GEN f) { return equali1(gel(f,1)); }
    4610             : 
    4611    15047085 : static GEN qfi_pow(void *E, GEN f, GEN n)
    4612    15047085 : { return E? nupow(f,n,(GEN)E): powgi(f,n); }
    4613    18997455 : static GEN qfi_comp(void *E, GEN f, GEN g)
    4614    18997455 : { return E? nucomp(f,g,(GEN)E): qficomp(f,g); }
    4615             : static const struct bb_group qfi_group={ qfi_comp,qfi_pow,NULL,hash_GEN,
    4616             :                                          gidentical,qfb_equal1,NULL};
    4617             : 
    4618             : GEN
    4619     2418537 : qfi_order(GEN q, GEN o)
    4620     2418537 : { return gen_order(q, o, NULL, &qfi_group); }
    4621             : 
    4622             : GEN
    4623           0 : qfi_log(GEN a, GEN g, GEN o)
    4624           0 : { return gen_PH_log(a, g, o, NULL, &qfi_group); }
    4625             : 
    4626             : GEN
    4627      521297 : qfi_Shanks(GEN a, GEN g, long n)
    4628             : {
    4629      521297 :   pari_sp av = avma;
    4630             :   GEN T, X;
    4631             :   long rt_n, c;
    4632             : 
    4633      521297 :   a = redimag(a);
    4634      521297 :   g = redimag(g);
    4635             : 
    4636      521297 :   rt_n = sqrt((double)n);
    4637      521297 :   c = n / rt_n;
    4638      521297 :   c = (c * rt_n < n + 1) ? c + 1 : c;
    4639             : 
    4640      521297 :   T = gen_Shanks_init(g, rt_n, NULL, &qfi_group);
    4641      521297 :   X = gen_Shanks(T, a, c, NULL, &qfi_group);
    4642             : 
    4643      521297 :   if (!X) { avma = av; return X; }
    4644      279692 :   return gerepileuptoint(av, X);
    4645             : }
    4646             : 
    4647             : GEN
    4648         140 : qfbclassno0(GEN x,long flag)
    4649             : {
    4650         140 :   switch(flag)
    4651             :   {
    4652         126 :     case 0: return map_proto_G(classno,x);
    4653          14 :     case 1: return map_proto_G(classno2,x);
    4654           0 :     default: pari_err_FLAG("qfbclassno");
    4655             :   }
    4656           0 :   return NULL; /* not reached */
    4657             : }
    4658             : 
    4659             : /* f^h = 1, return order(f). Set *pfao to its factorization */
    4660             : static GEN
    4661     2270664 : find_order(void *E, GEN f, GEN h, GEN *pfao)
    4662             : {
    4663     2270664 :   GEN v = gen_factored_order(f, h, E, &qfi_group);
    4664     2270664 :   *pfao = gel(v,2); return gel(v,1);
    4665             : }
    4666             : 
    4667             : static int
    4668        5550 : ok_q(GEN q, GEN h, GEN d2, long r2)
    4669             : {
    4670        5550 :   if (d2)
    4671             :   {
    4672           7 :     if (r2 <= 2 && !mpodd(q)) return 0;
    4673           7 :     return is_pm1(coprime_part(q,d2));
    4674             :   }
    4675             :   else
    4676             :   {
    4677        5543 :     if (r2 <= 1 && !mpodd(q)) return 0;
    4678        5543 :     return is_pm1(coprime_part(q,h));
    4679             :   }
    4680             : }
    4681             : 
    4682             : /* a,b given by their factorizations. Return factorization of lcm(a,b).
    4683             :  * Set A,B such that A*B = lcm(a, b), (A,B)=1, A|a, B|b */
    4684             : static GEN
    4685      295900 : split_lcm(GEN a, GEN Fa, GEN b, GEN Fb, GEN *pA, GEN *pB)
    4686             : {
    4687      295900 :   GEN P = ZV_union_shallow(gel(Fa,1), gel(Fb,1));
    4688      295900 :   GEN A = gen_1, B = gen_1;
    4689      295900 :   long i, l = lg(P);
    4690      295900 :   GEN E = cgetg(l, t_COL);
    4691      874865 :   for (i=1; i<l; i++)
    4692             :   {
    4693      578965 :     GEN p = gel(P,i);
    4694      578965 :     long va = Z_pval(a,p);
    4695      578965 :     long vb = Z_pval(b,p);
    4696      578965 :     if (va < vb)
    4697             :     {
    4698      297621 :       B = mulii(B,powiu(p,vb));
    4699      297621 :       gel(E,i) = utoi(vb);
    4700             :     }
    4701             :     else
    4702             :     {
    4703      281344 :       A = mulii(A,powiu(p,va));
    4704      281344 :       gel(E,i) = utoi(va);
    4705             :     }
    4706             :   }
    4707      295900 :   *pA = A;
    4708      295900 :   *pB = B; return mkmat2(P,E);
    4709             : }
    4710             : 
    4711             : /* g1 has order d1, f has order o, replace g1 by an element of order lcm(d1,o)*/
    4712             : static void
    4713      295900 : update_g1(GEN *pg1, GEN *pd1, GEN *pfad1, GEN f, GEN o, GEN fao)
    4714             : {
    4715      295900 :   GEN A,B, g1 = *pg1, d1 = *pd1;
    4716      295900 :   *pfad1 = split_lcm(d1,*pfad1, o,fao, &A,&B);
    4717      295900 :   *pg1 = gmul(powgi(g1, diviiexact(d1,A)),  powgi(f, diviiexact(o,B)));
    4718      295900 :   *pd1 = mulii(A,B); /* g1 has order d1 <- lcm(d1,o) */
    4719      295900 : }
    4720             : 
    4721             : /* Write x = Df^2, where D = fundamental discriminant,
    4722             :  * P^E = factorisation of conductor f, with E[i] >= 0 */
    4723             : static void
    4724     1720764 : corediscfact(GEN x, long xmod4, GEN *ptD, GEN *ptP, GEN *ptE)
    4725             : {
    4726     1720764 :   long s = signe(x), l, i;
    4727     1720764 :   GEN fa = absi_factor(x);
    4728     1720764 :   GEN d, P = gel(fa,1), E = gtovecsmall(gel(fa,2));
    4729             : 
    4730     1720764 :   l = lg(P); d = gen_1;
    4731     4508027 :   for (i=1; i<l; i++)
    4732             :   {
    4733     2787263 :     if (E[i] & 1) d = mulii(d, gel(P,i));
    4734     2787263 :     E[i] >>= 1;
    4735             :   }
    4736     1720764 :   if (!xmod4 && mod4(d) != ((s < 0)? 3: 1)) { d = shifti(d,2); E[1]--; }
    4737     1720764 :   *ptD = (s < 0)? negi(d): d;
    4738     1720764 :   *ptP = P;
    4739     1720764 :   *ptE = E;
    4740     1720764 : }
    4741             : 
    4742             : static GEN
    4743     1678710 : conductor_part(GEN x, long xmod4, GEN *ptD, GEN *ptreg)
    4744             : {
    4745     1678710 :   long l, i, s = signe(x);
    4746             :   GEN E, H, D, P, reg;
    4747             : 
    4748     1678710 :   corediscfact(x, xmod4, &D, &P, &E);
    4749     1678710 :   H = gen_1; l = lg(P);
    4750             :   /* f \prod_{p|f}  [ 1 - (D/p) p^-1 ] = \prod_{p^e||f} p^(e-1) [ p - (D/p) ] */
    4751     4360301 :   for (i=1; i<l; i++)
    4752             :   {
    4753     2681591 :     long e = E[i];
    4754     2681591 :     if (e)
    4755             :     {
    4756           7 :       GEN p = gel(P,i);
    4757           7 :       H = mulii(H, subis(p, kronecker(D,p)));
    4758           7 :       if (e >= 2) H = mulii(H, powiu(p,e-1));
    4759             :     }
    4760             :   }
    4761             : 
    4762             :   /* divide by [ O_K^* : O^* ] */
    4763     1678710 :   if (s < 0)
    4764             :   {
    4765     1678696 :     reg = NULL;
    4766     1678696 :     switch(itou_or_0(D))
    4767             :     {
    4768           0 :       case 4: H = shifti(H,-1); break;
    4769           0 :       case 3: H = divis(H,3); break;
    4770             :     }
    4771             :   } else {
    4772          14 :     reg = quadregulator(D,DEFAULTPREC);
    4773          14 :     if (!equalii(x,D))
    4774           0 :       H = divii(H, roundr(divrr(quadregulator(x,DEFAULTPREC), reg)));
    4775             :   }
    4776     1678710 :   if (ptreg) *ptreg = reg;
    4777     1678710 :   *ptD = D; return H;
    4778             : }
    4779             : 
    4780             : static long
    4781     1678689 : two_rank(GEN x)
    4782             : {
    4783     1678689 :   GEN p = gel(absi_factor(x),1);
    4784     1678689 :   long l = lg(p)-1;
    4785             : #if 0 /* positive disc not needed */
    4786             :   if (signe(x) > 0)
    4787             :   {
    4788             :     long i;
    4789             :     for (i=1; i<=l; i++)
    4790             :       if (mod4(gel(p,i)) == 3) { l--; break; }
    4791             :   }
    4792             : #endif
    4793     1678689 :   return l-1;
    4794             : }
    4795             : 
    4796             : static GEN
    4797    31893332 : sqr_primeform(GEN x, ulong p) { return redimag(qfisqr(primeform_u(x, p))); }
    4798             : /* return a set of forms hopefully generating Cl(K)^2; set L ~ L(chi_D,1) */
    4799             : static GEN
    4800     1678689 : get_forms(GEN D, GEN *pL)
    4801             : {
    4802     1678689 :   const long MAXFORM = 20;
    4803     1678689 :   GEN L, sqrtD = gsqrt(absi(D),DEFAULTPREC), forms = vectrunc_init(MAXFORM+1);
    4804     1678689 :   long s, nforms = 0;
    4805             :   ulong p;
    4806             :   forprime_t S;
    4807     1678689 :   L = mulrr(divrr(sqrtD,mppi(DEFAULTPREC)), dbltor(1.005));/*overshoot by 0.5%*/
    4808     1678689 :   s = itos_or_0( truncr(shiftr(sqrtr(sqrtD), 1)) );
    4809     1678689 :   if (!s) pari_err_OVERFLOW("classno [discriminant too large]");
    4810     1678689 :   if      (s < 10)   s = 200;
    4811     1554960 :   else if (s < 20)   s = 1000;
    4812        1477 :   else if (s < 5000) s = 5000;
    4813     1678689 :   u_forprime_init(&S, 2, s);
    4814   284027476 :   while ( (p = u_forprime_next(&S)) )
    4815             :   {
    4816   280670098 :     long d, k = kroiu(D,p);
    4817             :     pari_sp av2;
    4818   280670098 :     if (!k) continue;
    4819   278845499 :     if (k > 0)
    4820             :     {
    4821   139862100 :       if (++nforms < MAXFORM) vectrunc_append(forms, sqr_primeform(D,p));
    4822   139862100 :       d = p-1;
    4823             :     }
    4824             :     else
    4825   138983399 :       d = p+1;
    4826   278845499 :     av2 = avma; affrr(divru(mulur(p,L),d), L); avma = av2;
    4827             :   }
    4828     1678689 :   *pL = L; return forms;
    4829             : }
    4830             : 
    4831             : /* h ~ #G, return o = order of f, set fao = its factorization */
    4832             : static  GEN
    4833     1678738 : Shanks_order(void *E, GEN f, GEN h, GEN *pfao)
    4834             : {
    4835     1678738 :   long s = minss(itos(sqrti(h)), 10000);
    4836     1678738 :   GEN T = gen_Shanks_init(f, s, E, &qfi_group);
    4837     1678738 :   GEN v = gen_Shanks(T, ginv(f), ULONG_MAX, E, &qfi_group);
    4838     1678738 :   return find_order(E, f, addiu(v,1), pfao);
    4839             : }
    4840             : 
    4841             : /* if g = 1 in  G/<f> ? */
    4842             : static int
    4843         518 : equal1(void *E, GEN T, ulong N, GEN g)
    4844         518 : { return !!gen_Shanks(T, g, N, E, &qfi_group); }
    4845             : 
    4846             : /* Order of 'a' in G/<f>, T = gen_Shanks_init(f,n), order(f) < n*N
    4847             :  * FIXME: should be gen_order, but equal1 has the wrong prototype */
    4848             : static GEN
    4849         112 : relative_order(void *E, GEN a, GEN o, ulong N,  GEN T)
    4850             : {
    4851         112 :   pari_sp av = avma;
    4852             :   long i, l;
    4853             :   GEN m;
    4854             : 
    4855         112 :   m = dlog_get_ordfa(o);
    4856         112 :   if (!m) pari_err_TYPE("gen_order [missing order]",a);
    4857         112 :   o = gel(m,1);
    4858         112 :   m = gel(m,2); l = lgcols(m);
    4859         322 :   for (i = l-1; i; i--)
    4860             :   {
    4861         210 :     GEN t, y, p = gcoeff(m,i,1);
    4862         210 :     long j, e = itos(gcoeff(m,i,2));
    4863         210 :     if (l == 2) {
    4864          35 :       t = gen_1;
    4865          35 :       y = a;
    4866             :     } else {
    4867         175 :       t = diviiexact(o, powiu(p,e));
    4868         175 :       y = powgi(a, t);
    4869             :     }
    4870         210 :     if (equal1(E, T,N,y)) o = t;
    4871             :     else {
    4872         126 :       for (j = 1; j < e; j++)
    4873             :       {
    4874          28 :         y = powgi(y, p);
    4875          28 :         if (equal1(E, T,N,y)) break;
    4876             :       }
    4877         119 :       if (j < e) {
    4878          21 :         if (j > 1) p = powiu(p, j);
    4879          21 :         o = mulii(t, p);
    4880             :       }
    4881             :     }
    4882             :   }
    4883         112 :   return gerepilecopy(av, o);
    4884             : }
    4885             : 
    4886             : /* h(x) for x<0 using Baby Step/Giant Step.
    4887             :  * Assumes G is not too far from being cyclic.
    4888             :  *
    4889             :  * Compute G^2 instead of G so as to kill most of the non-cyclicity */
    4890             : GEN
    4891     1680638 : classno(GEN x)
    4892             : {
    4893     1680638 :   pari_sp av = avma;
    4894             :   long r2, k, s, i, l;
    4895             :   GEN forms, hin, Hf, D, g1, d1, d2, q, L, fad1, order_bound;
    4896             :   void *E;
    4897             : 
    4898     1680638 :   if (signe(x) >= 0) return classno2(x);
    4899             : 
    4900     1680631 :   check_quaddisc(x, &s, &k, "classno");
    4901     1680631 :   if (cmpiu(x,12) <= 0) return gen_1;
    4902             : 
    4903     1678689 :   Hf = conductor_part(x, k, &D, NULL);
    4904     1678689 :   if (cmpiu(D,12) <= 0) return gerepilecopy(av, Hf);
    4905     1678689 :   forms =  get_forms(D, &L);
    4906     1678689 :   r2 = two_rank(D);
    4907     1678689 :   hin = roundr(shiftr(L, -r2)); /* rough approximation for #G, G = Cl(K)^2 */
    4908             : 
    4909     1678689 :   l = lg(forms);
    4910     1678689 :   order_bound = const_vec(l-1, NULL);
    4911     1678689 :   E = expi(D) > 60? (void*)sqrtnint(shifti(absi(D),-2),4): NULL;
    4912     1678689 :   g1 = gel(forms,1);
    4913     1678689 :   gel(order_bound,1) = d1 = Shanks_order(E, g1, hin, &fad1);
    4914     1678689 :   q = diviiround(hin, d1); /* approximate order of G/<g1> */
    4915     1678689 :   d2 = NULL; /* not computed yet */
    4916     1678689 :   if (is_pm1(q)) goto END;
    4917      417553 :   for (i=2; i < l; i++)
    4918             :   {
    4919      411940 :     GEN o, fao, a, F, fd, f = gel(forms,i);
    4920      411940 :     fd = powgi(f, d1); if (is_pm1(gel(fd,1))) continue;
    4921      295900 :     F = powgi(fd, q);
    4922      295900 :     a = gel(F,1);
    4923      295900 :     o = is_pm1(a)? find_order(E, fd, q, &fao): Shanks_order(E, fd, q, &fao);
    4924             :     /* f^(d1 q) = 1 */
    4925      295900 :     fao = merge_factor(fad1,fao, (void*)&cmpii, &cmp_nodata);
    4926      295900 :     o = find_order(E, f, fao, &fao);
    4927      295900 :     gel(order_bound,i) = o;
    4928             :     /* o = order of f, fao = factor(o) */
    4929      295900 :     update_g1(&g1,&d1,&fad1, f,o,fao);
    4930      295900 :     q = diviiround(hin, d1);
    4931      295900 :     if (is_pm1(q)) goto END;
    4932             :   }
    4933             :   /* very probably d1 = expo(Cl^2(K)), q ~ #Cl^2(K) / d1 */
    4934        5613 :   if (expi(q) > 3)
    4935             :   { /* q large: compute d2, 2nd elt divisor */
    4936          70 :     ulong N, n = 2*itou(sqrti(d1));
    4937          70 :     GEN D = d1, T = gen_Shanks_init(g1, n, E, &qfi_group);
    4938          70 :     d2 = gen_1;
    4939          70 :     N = itou( gceil(gdivgs(d1,n)) ); /* order(g1) <= n*N */
    4940         287 :     for (i = 1; i < l; i++)
    4941             :     {
    4942         280 :       GEN d, f = gel(forms,i), B = gel(order_bound,i);
    4943         280 :       if (!B) B = find_order(E, f, fad1, /*junk*/&d);
    4944         280 :       f = powgi(f,d2);
    4945         280 :       if (equal1(E, T,N,f)) continue;
    4946         112 :       B = gdiv(B,d2); if (typ(B) == t_FRAC) B = gel(B,1);
    4947             :       /* f^B = 1 */
    4948         112 :       d = relative_order(E, f, B, N,T);
    4949         112 :       d2= mulii(d,d2);
    4950         112 :       D = mulii(d1,d2);
    4951         112 :       q = diviiround(hin,D);
    4952         112 :       if (is_pm1(q)) { d1 = D; goto END; }
    4953             :     }
    4954             :     /* very probably, d2 is the 2nd elementary divisor */
    4955           7 :     d1 = D; /* product of first two elt divisors */
    4956             :   }
    4957             :   /* impose q | d2^oo (d1^oo if d2 not computed), and compatible with known
    4958             :    * 2-rank */
    4959        5550 :   if (!ok_q(q,d1,d2,r2))
    4960             :   {
    4961           0 :     GEN q0 = q;
    4962             :     long d;
    4963           0 :     if (cmpii(mulii(q,d1), hin) < 0)
    4964             :     { /* try q = q0+1,-1,+2,-2 */
    4965           0 :       d = 1;
    4966           0 :       do { q = addis(q0,d); d = d>0? -d: 1-d; } while(!ok_q(q,d1,d2,r2));
    4967             :     }
    4968             :     else
    4969             :     { /* q0-1,+1,-2,+2  */
    4970           0 :       d = -1;
    4971           0 :       do { q = addis(q0,d); d = d<0? -d: -1-d; } while(!ok_q(q,d1,d2,r2));
    4972             :     }
    4973             :   }
    4974        5550 :   d1 = mulii(d1,q);
    4975             : 
    4976             : END:
    4977     1678689 :   return gerepileuptoint(av, shifti(mulii(d1,Hf), r2));
    4978             : }
    4979             : 
    4980             : /* use Euler products */
    4981             : GEN
    4982          21 : classno2(GEN x)
    4983             : {
    4984          21 :   pari_sp av = avma;
    4985          21 :   const long prec = DEFAULTPREC;
    4986             :   long n, i, r, s;
    4987             :   GEN p1, p2, S, p4, p5, p7, Hf, Pi, reg, logd, d, dr, D, half;
    4988             : 
    4989          21 :   check_quaddisc(x, &s, &r, "classno2");
    4990          21 :   if (s < 0 && cmpiu(x,12) <= 0) return gen_1;
    4991             : 
    4992          21 :   Hf = conductor_part(x, r, &D, &reg);
    4993          21 :   if (s < 0 && cmpiu(D,12) <= 0) return gerepilecopy(av, Hf); /* |D| < 12*/
    4994             : 
    4995          21 :   Pi = mppi(prec);
    4996          21 :   d = absi(D); dr = itor(d, prec);
    4997          21 :   logd = logr_abs(dr);
    4998          21 :   p1 = sqrtr(divrr(mulir(d,logd), gmul2n(Pi,1)));
    4999          21 :   if (s > 0)
    5000             :   {
    5001          14 :     GEN invlogd = invr(logd);
    5002          14 :     p2 = subsr(1, shiftr(mulrr(logr_abs(reg),invlogd),1));
    5003          14 :     if (cmprr(sqrr(p2), shiftr(invlogd,1)) >= 0) p1 = mulrr(p2,p1);
    5004             :   }
    5005          21 :   n = itos_or_0( mptrunc(p1) );
    5006          21 :   if (!n) pari_err_OVERFLOW("classno [discriminant too large]");
    5007             : 
    5008          21 :   p4 = divri(Pi,d);
    5009          21 :   p7 = invr(sqrtr_abs(Pi));
    5010          21 :   half = real2n(-1, prec);
    5011          21 :   if (s > 0)
    5012             :   { /* i = 1, shortcut */
    5013          14 :     p1 = sqrtr_abs(dr);
    5014          14 :     p5 = subsr(1, mulrr(p7,incgamc(half,p4,prec)));
    5015          14 :     S = addrr(mulrr(p1,p5), eint1(p4,prec));
    5016         546 :     for (i=2; i<=n; i++)
    5017             :     {
    5018         532 :       long k = kroiu(D,i); if (!k) continue;
    5019         434 :       p2 = mulir(sqru(i), p4);
    5020         434 :       p5 = subsr(1, mulrr(p7,incgamc(half,p2,prec)));
    5021         434 :       p5 = addrr(divru(mulrr(p1,p5),i), eint1(p2,prec));
    5022         434 :       S = (k>0)? addrr(S,p5): subrr(S,p5);
    5023             :     }
    5024          14 :     S = shiftr(divrr(S,reg),-1);
    5025             :   }
    5026             :   else
    5027             :   { /* i = 1, shortcut */
    5028           7 :     p1 = gdiv(sqrtr_abs(dr), Pi);
    5029           7 :     p5 = subsr(1, mulrr(p7,incgamc(half,p4,prec)));
    5030           7 :     S = addrr(p5, divrr(p1, mpexp(p4)));
    5031         952 :     for (i=2; i<=n; i++)
    5032             :     {
    5033         945 :       long k = kroiu(D,i); if (!k) continue;
    5034         945 :       p2 = mulir(sqru(i), p4);
    5035         945 :       p5 = subsr(1, mulrr(p7,incgamc(half,p2,prec)));
    5036         945 :       p5 = addrr(p5, divrr(p1, mulur(i, mpexp(p2))));
    5037         945 :       S = (k>0)? addrr(S,p5): subrr(S,p5);
    5038             :     }
    5039             :   }
    5040          21 :   return gerepileuptoint(av, mulii(Hf, roundr(S)));
    5041             : }
    5042             : 
    5043             : static GEN
    5044        5304 : hclassno2(GEN x)
    5045             : {
    5046             :   long i, l, s, xmod4;
    5047             :   GEN Q, H, D, P, E;
    5048             : 
    5049        5304 :   x = negi(x);
    5050        5304 :   check_quaddisc(x, &s, &xmod4, "hclassno");
    5051        5304 :   corediscfact(x, xmod4, &D, &P, &E);
    5052             : 
    5053        5304 :   Q = quadclassunit0(D, 0, NULL, 0);
    5054        5304 :   H = gel(Q,1); l = lg(P);
    5055             : 
    5056             :   /* H \prod_{p^e||f}  (1 + (p^e-1)/(p-1))[ p - (D/p) ] */
    5057       24043 :   for (i=1; i<l; i++)
    5058             :   {
    5059       18739 :     long e = E[i];
    5060       18739 :     if (e)
    5061             :     {
    5062        4198 :       GEN p = gel(P,i), t = subis(p, kronecker(D,p));
    5063        4198 :       if (e > 1) t = mulii(t, diviiexact(subis(powiu(p,e), 1), subis(p,1)));
    5064        4198 :       H = mulii(H, addsi(1, t));
    5065             :     }
    5066             :   }
    5067        5304 :   switch( itou_or_0(D) )
    5068             :   {
    5069           0 :     case 3: H = gdivgs(H, 3); break;
    5070           0 :     case 4: H = gdivgs(H, 2); break;
    5071             :   }
    5072        5304 :   return H;
    5073             : }
    5074             : 
    5075             : GEN
    5076       76111 : hclassno(GEN x)
    5077             : {
    5078             :   ulong a, b, b2, d, h;
    5079             :   long s;
    5080             :   int f;
    5081             : 
    5082       76111 :   if (typ(x) != t_INT) pari_err_TYPE("hclassno",x);
    5083       76111 :   s = signe(x);
    5084       76111 :   if (s < 0) return gen_0;
    5085       76111 :   if (!s) return gdivgs(gen_1, -12);
    5086             : 
    5087       76111 :   a = mod4(x); if (a == 1 || a == 2) return gen_0;
    5088             : 
    5089       76111 :   d = itou_or_0(x);
    5090       76111 :   if (!d || d > 500000) return hclassno2(x);
    5091             : 
    5092       70807 :   h = 0; b = d&1; b2 = (1+d)>>2; f=0;
    5093       70807 :   if (!b)
    5094             :   {
    5095     2074994 :     for (a=1; a*a<b2; a++)
    5096     2020913 :       if (b2%a == 0) h++;
    5097       54081 :     f = (a*a==b2); b=2; b2=(4+d)>>2;
    5098             :   }
    5099     2092598 :   while (b2*3 < d)
    5100             :   {
    5101     1950984 :     if (b2%b == 0) h++;
    5102   131085034 :     for (a=b+1; a*a < b2; a++)
    5103   129134050 :       if (b2%a == 0) h += 2;
    5104     1950984 :     if (a*a == b2) h++;
    5105     1950984 :     b += 2; b2 = (b*b+d)>>2;
    5106             :   }
    5107       70807 :   if (b2*3 == d) retmkfrac(utoipos(3*h+1), utoipos(3));
    5108       70702 :   if (f) retmkfrac(utoipos(2*h+1), gen_2);
    5109       66131 :   return utoipos(h);
    5110             : }
    5111             : 
    5112             : /******************************************************************/
    5113             : /*                                                                */
    5114             : /*                 RAMANUJAN's TAU FUNCTION                       */
    5115             : /*                                                                */
    5116             : /******************************************************************/
    5117             : 
    5118             : /* h(D) / (w(D)/2), D fundamental */
    5119             : static GEN
    5120       36750 : myh(GEN D)
    5121             : {
    5122             :   GEN Q;
    5123       36750 :   if (equalis(D, -3)) return mkfrac(gen_1, utoipos(3));
    5124       32690 :   if (equalis(D, -4)) return ghalf;
    5125       28161 :   Q = quadclassunit0(D, 0, NULL, 0); return gel(Q,1);
    5126             : }
    5127             : 
    5128             : /* 1 + q + ... + q^v */
    5129             : static GEN
    5130        6083 : geomsum(GEN q, long v)
    5131             : {
    5132             :   GEN u;
    5133        6083 :   if (v == 0) return gen_1;
    5134         273 :   u = addui(1, q);
    5135         273 :   for (; v > 1; v--) u = addui(1, mulii(q, u));
    5136         273 :   return u;
    5137             : }
    5138             : 
    5139             : /* 4|N, not fundamental at 2; Hurwitz class number in level 2,
    5140             :  * equal to H(N)+2H(N/4), H=qfbhclassno */
    5141             : static GEN
    5142       36750 : Hspec(GEN N)
    5143             : {
    5144             :   GEN D0, P, E, H, s;
    5145             :   long j, lP;
    5146             : 
    5147       36750 :   corediscfact(negi(N), 0, &D0, &P, &E);
    5148       36750 :   H = myh(D0);
    5149       36750 :   lP = lg(P); /* E[1] > 0 since N not fundamental at */
    5150       36750 :   s = addui(3, muliu(subiu(int2n(E[1]+1), 3), 2 - kroiu(D0,2)));
    5151       86933 :   for (j = 2; j < lP; j++)
    5152             :   {
    5153       50183 :     long e = E[j];
    5154       50183 :     if (e)
    5155             :     {
    5156        6083 :       ulong p = itou(gel(P,j));
    5157        6083 :       GEN q = geomsum(utoipos(p), e-1);
    5158        6083 :       s = mulii(s, addui(1, muliu(q, p - kroiu(D0,p))));
    5159             :     }
    5160             :   }
    5161       36750 :   return gmul(H,s);
    5162             : }
    5163             : 
    5164             : /* Ramanujan tau function for p prime */
    5165             : static GEN
    5166       14903 : tauprime(GEN p)
    5167             : {
    5168       14903 :   pari_sp av = avma, av2;
    5169             :   GEN s, p_27, p_9, T;
    5170             :   ulong lim, t, tin;
    5171             : 
    5172       14903 :   if (equaliu(p, 2)) return utoineg(24);
    5173             :   /* p > 2 */
    5174       11396 :   p_27 = mulsi(7, sqri(p));
    5175       11396 :   p_9 = mulsi(9, p);
    5176       11396 :   av2 = avma;
    5177       11396 :   lim = itou(sqrtint(p));
    5178       11396 :   tin = mod4(p) == 3? 1: 0;
    5179       11396 :   s = gen_0;
    5180       87178 :   for (t = 1; t <= lim; ++t)
    5181             :   {
    5182       75782 :     GEN p2, t2 = sqru(t), t3 = shifti(subii(p, t2), 2); /* 4(p-t^2) */
    5183             :     /* t mod 2 != tin <=> t3 not fundamental at 2 */
    5184       75782 :     p2 = ((t&1) == tin) ? hclassno(t3) : Hspec(t3);
    5185       75782 :     s = gadd(s, gmul(mulii(powiu(t2,3), addii(p_27, mulii(t2, subii(mulsi(4, t2), p_9)))), p2));
    5186       75782 :     if (!(t & 255)) s = gerepileupto(av2, s);
    5187             :   }
    5188             :   /* 28p^3 - 28p^2 - 90p - 35 */
    5189       11396 :   T = addsi(-35, mulii(p, addsi(-90, mulii(p, addsi(-28, mulsi(28, p))))));
    5190       11396 :   return gerepileupto(av, subii(mulii(powiu(p,3),T), addsi(1, gmulsg(128, s))));
    5191             : }
    5192             : 
    5193             : /* Ramanujan tau function, return 0 for <= 0 */
    5194             : GEN
    5195        7035 : ramanujantau(GEN n)
    5196             : {
    5197        7035 :   pari_sp ltop = avma;
    5198             :   GEN T, F, P, E;
    5199             :   long j, lP;
    5200             : 
    5201        7035 :   if (!(F = check_arith_all(n,"ramanujantau")))
    5202             :   {
    5203        7014 :     if (signe(n) <= 0) return gen_0;
    5204        7007 :     F = Z_factor(n);
    5205             :   }
    5206             :   else
    5207             :   {
    5208          21 :     P = gel(F,1);
    5209          21 :     if (lg(P) == 1 || signe(gel(P,1)) <= 0) return gen_0;
    5210             :   }
    5211             : 
    5212        7014 :   P = gel(F,1);
    5213        7014 :   E = gel(F,2); lP = lg(P);
    5214        7014 :   T = gen_1;
    5215       21917 :   for (j = 1; j < lP; j++)
    5216             :   {
    5217       14903 :     GEN p = gel(P,j), tp = tauprime(p), t1 = tp, t0 = gen_1;
    5218       14903 :     long k, e = itou(gel(E,j));
    5219       20160 :     for (k = 1; k < e; k++)
    5220             :     {
    5221        5257 :       GEN t2 = subii(mulii(tp, t1), mulii(powiu(p, 11), t0));
    5222        5257 :       t0 = t1; t1 = t2;
    5223             :     }
    5224       14903 :     T = mulii(T, t1);
    5225             :   }
    5226        7014 :   return gerepileuptoint(ltop, T);
    5227             : }

Generated by: LCOV version 1.11