Code coverage tests

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

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

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

LCOV - code coverage report
Current view: top level - basemath - arith1.c (source / functions) Hit Total Coverage
Test: PARI/GP v2.16.2 lcov report (development 29115-f22e516b23) Lines: 2099 2310 90.9 %
Date: 2024-03-18 08:03:28 Functions: 211 225 93.8 %
Legend: Lines: hit not hit

          Line data    Source code
       1             : /* Copyright (C) 2000  The PARI group.
       2             : 
       3             : This file is part of the PARI/GP package.
       4             : 
       5             : PARI/GP is free software; you can redistribute it and/or modify it under the
       6             : terms of the GNU General Public License as published by the Free Software
       7             : Foundation; either version 2 of the License, or (at your option) any later
       8             : version. It is distributed in the hope that it will be useful, but WITHOUT
       9             : ANY WARRANTY WHATSOEVER.
      10             : 
      11             : Check the License for details. You should have received a copy of it, along
      12             : with the package; see the file 'COPYING'. If not, write to the Free Software
      13             : Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. */
      14             : 
      15             : /*********************************************************************/
      16             : /**                     ARITHMETIC FUNCTIONS                        **/
      17             : /**                         (first part)                            **/
      18             : /*********************************************************************/
      19             : #include "pari.h"
      20             : #include "paripriv.h"
      21             : 
      22             : #define DEBUGLEVEL DEBUGLEVEL_arith
      23             : 
      24             : /******************************************************************/
      25             : /*                 GENERATOR of (Z/mZ)*                           */
      26             : /******************************************************************/
      27             : static GEN
      28        1812 : remove2(GEN q) { long v = vali(q); return v? shifti(q, -v): q; }
      29             : static ulong
      30      465082 : u_remove2(ulong q) { return q >> vals(q); }
      31             : GEN
      32        1812 : odd_prime_divisors(GEN q) { return gel(Z_factor(remove2(q)), 1); }
      33             : static GEN
      34      465082 : u_odd_prime_divisors(ulong q) { return gel(factoru(u_remove2(q)), 1); }
      35             : /* p odd prime, q=(p-1)/2; L0 list of (some) divisors of q = (p-1)/2 or NULL
      36             :  * (all prime divisors of q); return the q/l, l in L0 */
      37             : static GEN
      38        4909 : is_gener_expo(GEN p, GEN L0)
      39             : {
      40        4909 :   GEN L, q = shifti(p,-1);
      41             :   long i, l;
      42        4909 :   if (L0) {
      43        3134 :     l = lg(L0);
      44        3134 :     L = cgetg(l, t_VEC);
      45             :   } else {
      46        1775 :     L0 = L = odd_prime_divisors(q);
      47        1775 :     l = lg(L);
      48             :   }
      49       14199 :   for (i=1; i<l; i++) gel(L,i) = diviiexact(q, gel(L0,i));
      50        4909 :   return L;
      51             : }
      52             : static GEN
      53      532882 : u_is_gener_expo(ulong p, GEN L0)
      54             : {
      55      532882 :   const ulong q = p >> 1;
      56             :   long i;
      57             :   GEN L;
      58      532882 :   if (!L0) L0 = u_odd_prime_divisors(q);
      59      532882 :   L = cgetg_copy(L0,&i);
      60     1155340 :   while (--i) L[i] = q / uel(L0,i);
      61      532881 :   return L;
      62             : }
      63             : 
      64             : int
      65     1629741 : is_gener_Fl(ulong x, ulong p, ulong p_1, GEN L)
      66             : {
      67             :   long i;
      68     1629741 :   if (krouu(x, p) >= 0) return 0;
      69     1380307 :   for (i=lg(L)-1; i; i--)
      70             :   {
      71      842641 :     ulong t = Fl_powu(x, uel(L,i), p);
      72      842643 :     if (t == p_1 || t == 1) return 0;
      73             :   }
      74      537666 :   return 1;
      75             : }
      76             : /* assume p prime */
      77             : ulong
      78     1054452 : pgener_Fl_local(ulong p, GEN L0)
      79             : {
      80     1054452 :   const pari_sp av = avma;
      81     1054452 :   const ulong p_1 = p-1;
      82             :   long x;
      83             :   GEN L;
      84     1054452 :   if (p <= 19) switch(p)
      85             :   { /* quick trivial cases */
      86          28 :     case 2:  return 1;
      87      116204 :     case 7:
      88      116204 :     case 17: return 3;
      89      405373 :     default: return 2;
      90             :   }
      91      532847 :   L = u_is_gener_expo(p,L0);
      92     1622070 :   for (x = 2;; x++)
      93     1622070 :     if (is_gener_Fl(x,p,p_1,L)) return gc_ulong(av, x);
      94             : }
      95             : ulong
      96      579470 : pgener_Fl(ulong p) { return pgener_Fl_local(p, NULL); }
      97             : 
      98             : /* L[i] = set of (p-1)/2l, l ODD prime divisor of p-1 (l=2 can be included,
      99             :  * but wasteful) */
     100             : int
     101       13899 : is_gener_Fp(GEN x, GEN p, GEN p_1, GEN L)
     102             : {
     103       13899 :   long i, t = lgefint(x)==3? kroui(x[2], p): kronecker(x, p);
     104       13899 :   if (t >= 0) return 0;
     105       21929 :   for (i = lg(L)-1; i; i--)
     106             :   {
     107       14400 :     GEN t = Fp_pow(x, gel(L,i), p);
     108       14400 :     if (equalii(t, p_1) || equali1(t)) return 0;
     109             :   }
     110        7529 :   return 1;
     111             : }
     112             : 
     113             : /* assume p prime, return a generator of all L[i]-Sylows in F_p^*. */
     114             : GEN
     115      357862 : pgener_Fp_local(GEN p, GEN L0)
     116             : {
     117      357862 :   pari_sp av0 = avma;
     118             :   GEN x, p_1, L;
     119      357862 :   if (lgefint(p) == 3)
     120             :   {
     121             :     ulong z;
     122      352958 :     if (p[2] == 2) return gen_1;
     123      258032 :     if (L0) L0 = ZV_to_nv(L0);
     124      258035 :     z = pgener_Fl_local(uel(p,2), L0);
     125      258056 :     return gc_utoipos(av0, z);
     126             :   }
     127        4904 :   p_1 = subiu(p,1); L = is_gener_expo(p, L0);
     128        4904 :   x = utoipos(2);
     129        9927 :   for (;; x[2]++) { if (is_gener_Fp(x, p, p_1, L)) break; }
     130        4904 :   return gc_utoipos(av0, uel(x,2));
     131             : }
     132             : 
     133             : GEN
     134       44135 : pgener_Fp(GEN p) { return pgener_Fp_local(p, NULL); }
     135             : 
     136             : ulong
     137      203414 : pgener_Zl(ulong p)
     138             : {
     139      203414 :   if (p == 2) pari_err_DOMAIN("pgener_Zl","p","=",gen_2,gen_2);
     140             :   /* only p < 2^32 such that znprimroot(p) != znprimroot(p^2) */
     141      203414 :   if (p == 40487) return 10;
     142             : #ifndef LONG_IS_64BIT
     143       29480 :   return pgener_Fl(p);
     144             : #else
     145      173934 :   if (p < (1UL<<32)) return pgener_Fl(p);
     146             :   else
     147             :   {
     148          30 :     const pari_sp av = avma;
     149          30 :     const ulong p_1 = p-1;
     150             :     long x ;
     151          30 :     GEN p2 = sqru(p), L = u_is_gener_expo(p, NULL);
     152         102 :     for (x=2;;x++)
     153         102 :       if (is_gener_Fl(x,p,p_1,L) && !is_pm1(Fp_powu(utoipos(x),p_1,p2)))
     154          30 :         return gc_ulong(av, x);
     155             :   }
     156             : #endif
     157             : }
     158             : 
     159             : /* p prime. Return a primitive root modulo p^e, e > 1 */
     160             : GEN
     161      168665 : pgener_Zp(GEN p)
     162             : {
     163      168665 :   if (lgefint(p) == 3) return utoipos(pgener_Zl(p[2]));
     164             :   else
     165             :   {
     166           5 :     const pari_sp av = avma;
     167           5 :     GEN p_1 = subiu(p,1), p2 = sqri(p), L = is_gener_expo(p,NULL);
     168           5 :     GEN x = utoipos(2);
     169          12 :     for (;; x[2]++)
     170          17 :       if (is_gener_Fp(x,p,p_1,L) && !equali1(Fp_pow(x,p_1,p2))) break;
     171           5 :     return gc_utoipos(av, uel(x,2));
     172             :   }
     173             : }
     174             : 
     175             : static GEN
     176         259 : gener_Zp(GEN q, GEN F)
     177             : {
     178         259 :   GEN p = NULL;
     179         259 :   long e = 0;
     180         259 :   if (F)
     181             :   {
     182          14 :     GEN P = gel(F,1), E = gel(F,2);
     183          14 :     long i, l = lg(P);
     184          42 :     for (i = 1; i < l; i++)
     185             :     {
     186          28 :       p = gel(P,i);
     187          28 :       if (absequaliu(p, 2)) continue;
     188          14 :       if (i < l-1) pari_err_DOMAIN("znprimroot", "n","=",F,F);
     189          14 :       e = itos(gel(E,i));
     190             :     }
     191          14 :     if (!p) pari_err_DOMAIN("znprimroot", "n","=",F,F);
     192             :   }
     193             :   else
     194         245 :     e = Z_isanypower(q, &p);
     195         259 :   if (!BPSW_psp(e? p: q)) pari_err_DOMAIN("znprimroot", "n","=", q,q);
     196         245 :   return e > 1? pgener_Zp(p): pgener_Fp(q);
     197             : }
     198             : 
     199             : GEN
     200         329 : znprimroot(GEN N)
     201             : {
     202         329 :   pari_sp av = avma;
     203             :   GEN x, n, F;
     204             : 
     205         329 :   if ((F = check_arith_non0(N,"znprimroot")))
     206             :   {
     207          14 :     F = clean_Z_factor(F);
     208          14 :     N = typ(N) == t_VEC? gel(N,1): factorback(F);
     209             :   }
     210         322 :   N = absi_shallow(N);
     211         322 :   if (abscmpiu(N, 4) <= 0) { set_avma(av); return mkintmodu(N[2]-1,N[2]); }
     212         273 :   switch(mod4(N))
     213             :   {
     214          14 :     case 0: /* N = 0 mod 4 */
     215          14 :       pari_err_DOMAIN("znprimroot", "n","=",N,N);
     216           0 :       x = NULL; break;
     217          28 :     case 2: /* N = 2 mod 4 */
     218          28 :       n = shifti(N,-1); /* becomes odd */
     219          28 :       x = gener_Zp(n,F); if (!mod2(x)) x = addii(x,n);
     220          21 :       break;
     221         231 :     default: /* N odd */
     222         231 :       x = gener_Zp(N,F);
     223         224 :       break;
     224             :   }
     225         245 :   return gerepilecopy(av, mkintmod(x, N));
     226             : }
     227             : 
     228             : /* n | (p-1), returns a primitive n-th root of 1 in F_p^* */
     229             : GEN
     230           0 : rootsof1_Fp(GEN n, GEN p)
     231             : {
     232           0 :   pari_sp av = avma;
     233           0 :   GEN L = odd_prime_divisors(n); /* 2 implicit in pgener_Fp_local */
     234           0 :   GEN z = pgener_Fp_local(p, L);
     235           0 :   z = Fp_pow(z, diviiexact(subiu(p,1), n), p); /* prim. n-th root of 1 */
     236           0 :   return gerepileuptoint(av, z);
     237             : }
     238             : 
     239             : GEN
     240        3033 : rootsof1u_Fp(ulong n, GEN p)
     241             : {
     242        3033 :   pari_sp av = avma;
     243        3033 :   GEN z, L = u_odd_prime_divisors(n); /* 2 implicit in pgener_Fp_local */
     244        3033 :   z = pgener_Fp_local(p, Flv_to_ZV(L));
     245        3033 :   z = Fp_pow(z, diviuexact(subiu(p,1), n), p); /* prim. n-th root of 1 */
     246        3033 :   return gerepileuptoint(av, z);
     247             : }
     248             : 
     249             : ulong
     250      214770 : rootsof1_Fl(ulong n, ulong p)
     251             : {
     252      214770 :   pari_sp av = avma;
     253      214770 :   GEN L = u_odd_prime_divisors(n); /* 2 implicit in pgener_Fl_local */
     254      214770 :   ulong z = pgener_Fl_local(p, L);
     255      214770 :   z = Fl_powu(z, (p-1) / n, p); /* prim. n-th root of 1 */
     256      214770 :   return gc_ulong(av,z);
     257             : }
     258             : 
     259             : /*********************************************************************/
     260             : /**                     INVERSE TOTIENT FUNCTION                    **/
     261             : /*********************************************************************/
     262             : /* N t_INT, L a ZV containing all prime divisors of N, and possibly other
     263             :  * primes. Return factor(N) */
     264             : GEN
     265      350651 : Z_factor_listP(GEN N, GEN L)
     266             : {
     267      350651 :   long i, k, l = lg(L);
     268      350651 :   GEN P = cgetg(l, t_COL), E = cgetg(l, t_COL);
     269     1346688 :   for (i = k = 1; i < l; i++)
     270             :   {
     271      996037 :     GEN p = gel(L,i);
     272      996037 :     long v = Z_pvalrem(N, p, &N);
     273      996037 :     if (v)
     274             :     {
     275      792176 :       gel(P,k) = p;
     276      792176 :       gel(E,k) = utoipos(v);
     277      792176 :       k++;
     278             :     }
     279             :   }
     280      350651 :   setlg(P, k);
     281      350651 :   setlg(E, k); return mkmat2(P,E);
     282             : }
     283             : 
     284             : /* look for x such that phi(x) = n, p | x => p > m (if m = NULL: no condition).
     285             :  * L is a list of primes containing all prime divisors of n. */
     286             : static long
     287      621565 : istotient_i(GEN n, GEN m, GEN L, GEN *px)
     288             : {
     289      621565 :   pari_sp av = avma, av2;
     290             :   GEN k, D;
     291             :   long i, v;
     292      621565 :   if (m && mod2(n))
     293             :   {
     294      270914 :     if (!equali1(n)) return 0;
     295       69986 :     if (px) *px = gen_1;
     296       69986 :     return 1;
     297             :   }
     298      350651 :   D = divisors(Z_factor_listP(shifti(n, -1), L));
     299             :   /* loop through primes p > m, d = p-1 | n */
     300      350651 :   av2 = avma;
     301      350651 :   if (!m)
     302             :   { /* special case p = 2, d = 1 */
     303       69986 :     k = n;
     304       69986 :     for (v = 1;; v++) {
     305       69986 :       if (istotient_i(k, gen_2, L, px)) {
     306       69986 :         if (px) *px = shifti(*px, v);
     307       69986 :         return 1;
     308             :       }
     309           0 :       if (mod2(k)) break;
     310           0 :       k = shifti(k,-1);
     311             :     }
     312           0 :     set_avma(av2);
     313             :   }
     314     1099462 :   for (i = 1; i < lg(D); ++i)
     315             :   {
     316     1001588 :     GEN p, d = shifti(gel(D, i), 1); /* even divisors of n */
     317     1001588 :     if (m && cmpii(d, m) < 0) continue;
     318      677782 :     p = addiu(d, 1);
     319      677782 :     if (!isprime(p)) continue;
     320      442064 :     k = diviiexact(n, d);
     321      481593 :     for (v = 1;; v++) {
     322             :       GEN r;
     323      481593 :       if (istotient_i(k, p, L, px)) {
     324      182791 :         if (px) *px = mulii(*px, powiu(p, v));
     325      182791 :         return 1;
     326             :       }
     327      298802 :       k = dvmdii(k, p, &r);
     328      298802 :       if (r != gen_0) break;
     329             :     }
     330      259273 :     set_avma(av2);
     331             :   }
     332       97874 :   return gc_long(av,0);
     333             : }
     334             : 
     335             : /* find x such that phi(x) = n */
     336             : long
     337       70000 : istotient(GEN n, GEN *px)
     338             : {
     339       70000 :   pari_sp av = avma;
     340       70000 :   if (typ(n) != t_INT) pari_err_TYPE("istotient", n);
     341       70000 :   if (signe(n) < 1) return 0;
     342       70000 :   if (mod2(n))
     343             :   {
     344          14 :     if (!equali1(n)) return 0;
     345          14 :     if (px) *px = gen_1;
     346          14 :     return 1;
     347             :   }
     348       69986 :   if (istotient_i(n, NULL, gel(Z_factor(n), 1), px))
     349             :   {
     350       69986 :     if (!px) set_avma(av);
     351             :     else
     352       69986 :       *px = gerepileuptoint(av, *px);
     353       69986 :     return 1;
     354             :   }
     355           0 :   return gc_long(av,0);
     356             : }
     357             : 
     358             : /*********************************************************************/
     359             : /**                        KRONECKER SYMBOL                         **/
     360             : /*********************************************************************/
     361             : /* t = 3,5 mod 8 ?  (= 2 not a square mod t) */
     362             : static int
     363   317684136 : ome(long t)
     364             : {
     365   317684136 :   switch(t & 7)
     366             :   {
     367   180017509 :     case 3:
     368   180017509 :     case 5: return 1;
     369   137666627 :     default: return 0;
     370             :   }
     371             : }
     372             : /* t a t_INT, is t = 3,5 mod 8 ? */
     373             : static int
     374     5582965 : gome(GEN t)
     375     5582965 : { return signe(t)? ome( mod2BIL(t) ): 0; }
     376             : 
     377             : /* assume y odd, return kronecker(x,y) * s */
     378             : static long
     379   223992413 : krouu_s(ulong x, ulong y, long s)
     380             : {
     381   223992413 :   ulong x1 = x, y1 = y, z;
     382  1018734920 :   while (x1)
     383             :   {
     384   794780580 :     long r = vals(x1);
     385   794830988 :     if (r)
     386             :     {
     387   423041591 :       if (odd(r) && ome(y1)) s = -s;
     388   422953110 :       x1 >>= r;
     389             :     }
     390   794742507 :     if (x1 & y1 & 2) s = -s;
     391   794742507 :     z = y1 % x1; y1 = x1; x1 = z;
     392             :   }
     393   223954340 :   return (y1 == 1)? s: 0;
     394             : }
     395             : 
     396             : long
     397    11938845 : kronecker(GEN x, GEN y)
     398             : {
     399    11938845 :   pari_sp av = avma;
     400    11938845 :   long s = 1, r;
     401             :   ulong xu;
     402             : 
     403    11938845 :   if (typ(x) != t_INT) pari_err_TYPE("kronecker",x);
     404    11938845 :   if (typ(y) != t_INT) pari_err_TYPE("kronecker",y);
     405    11938845 :   switch (signe(y))
     406             :   {
     407          63 :     case -1: y = negi(y); if (signe(x) < 0) s = -1; break;
     408         133 :     case 0: return is_pm1(x);
     409             :   }
     410    11938712 :   r = vali(y);
     411    11938709 :   if (r)
     412             :   {
     413     1348897 :     if (!mpodd(x)) return gc_long(av,0);
     414      321710 :     if (odd(r) && gome(x)) s = -s;
     415      321710 :     y = shifti(y,-r);
     416             :   }
     417    10911522 :   x = modii(x,y);
     418    13294729 :   while (lgefint(x) > 3) /* x < y */
     419             :   {
     420             :     GEN z;
     421     2383298 :     r = vali(x);
     422     2382995 :     if (r)
     423             :     {
     424     1302127 :       if (odd(r) && gome(y)) s = -s;
     425     1302056 :       x = shifti(x,-r);
     426             :     }
     427             :     /* x=3 mod 4 && y=3 mod 4 ? (both are odd here) */
     428     2382546 :     if (mod2BIL(x) & mod2BIL(y) & 2) s = -s;
     429     2381957 :     z = remii(y,x); y = x; x = z;
     430     2383390 :     if (gc_needed(av,2))
     431             :     {
     432           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"kronecker");
     433           0 :       gerepileall(av, 2, &x, &y);
     434             :     }
     435             :   }
     436    10911431 :   xu = itou(x);
     437    10911427 :   if (!xu) return is_pm1(y)? s: 0;
     438    10814235 :   r = vals(xu);
     439    10814234 :   if (r)
     440             :   {
     441     5744456 :     if (odd(r) && gome(y)) s = -s;
     442     5744456 :     xu >>= r;
     443             :   }
     444             :   /* x=3 mod 4 && y=3 mod 4 ? (both are odd here) */
     445    10814234 :   if (xu & mod2BIL(y) & 2) s = -s;
     446    10814237 :   return gc_long(av, krouu_s(umodiu(y,xu), xu, s));
     447             : }
     448             : 
     449             : long
     450       31157 : krois(GEN x, long y)
     451             : {
     452             :   ulong yu;
     453       31157 :   long s = 1;
     454             : 
     455       31157 :   if (y <= 0)
     456             :   {
     457          28 :     if (y == 0) return is_pm1(x);
     458           0 :     yu = (ulong)-y; if (signe(x) < 0) s = -1;
     459             :   }
     460             :   else
     461       31129 :     yu = (ulong)y;
     462       31129 :   if (!odd(yu))
     463             :   {
     464             :     long r;
     465       14819 :     if (!mpodd(x)) return 0;
     466       10913 :     r = vals(yu); yu >>= r;
     467       10913 :     if (odd(r) && gome(x)) s = -s;
     468             :   }
     469       27223 :   return krouu_s(umodiu(x, yu), yu, s);
     470             : }
     471             : /* assume y != 0 */
     472             : long
     473    27324259 : kroiu(GEN x, ulong y)
     474             : {
     475             :   long r;
     476    27324259 :   if (odd(y)) return krouu_s(umodiu(x,y), y, 1);
     477      299052 :   if (!mpodd(x)) return 0;
     478      204479 :   r = vals(y); y >>= r;
     479      204481 :   return krouu_s(umodiu(x,y), y, (odd(r) && gome(x))? -1: 1);
     480             : }
     481             : 
     482             : /* assume y > 0, odd, return s * kronecker(x,y) */
     483             : static long
     484      179211 : krouodd(ulong x, GEN y, long s)
     485             : {
     486             :   long r;
     487      179211 :   if (lgefint(y) == 3) return krouu_s(x, y[2], s);
     488       29273 :   if (!x) return 0; /* y != 1 */
     489       29273 :   r = vals(x);
     490       29273 :   if (r)
     491             :   {
     492       15242 :     if (odd(r) && gome(y)) s = -s;
     493       15242 :     x >>= r;
     494             :   }
     495             :   /* x=3 mod 4 && y=3 mod 4 ? (both are odd here) */
     496       29273 :   if (x & mod2BIL(y) & 2) s = -s;
     497       29273 :   return krouu_s(umodiu(y,x), x, s);
     498             : }
     499             : 
     500             : long
     501      144257 : krosi(long x, GEN y)
     502             : {
     503      144257 :   const pari_sp av = avma;
     504      144257 :   long s = 1, r;
     505      144257 :   switch (signe(y))
     506             :   {
     507           0 :     case -1: y = negi(y); if (x < 0) s = -1; break;
     508           0 :     case 0: return (x==1 || x==-1);
     509             :   }
     510      144257 :   r = vali(y);
     511      144257 :   if (r)
     512             :   {
     513       16884 :     if (!odd(x)) return gc_long(av,0);
     514       16884 :     if (odd(r) && ome(x)) s = -s;
     515       16884 :     y = shifti(y,-r);
     516             :   }
     517      144257 :   if (x < 0) { x = -x; if (mod4(y) == 3) s = -s; }
     518      144257 :   return gc_long(av, krouodd((ulong)x, y, s));
     519             : }
     520             : 
     521             : long
     522       34954 : kroui(ulong x, GEN y)
     523             : {
     524       34954 :   const pari_sp av = avma;
     525       34954 :   long s = 1, r;
     526       34954 :   switch (signe(y))
     527             :   {
     528           0 :     case -1: y = negi(y); break;
     529           0 :     case 0: return x==1UL;
     530             :   }
     531       34954 :   r = vali(y);
     532       34954 :   if (r)
     533             :   {
     534           0 :     if (!odd(x)) return gc_long(av,0);
     535           0 :     if (odd(r) && ome(x)) s = -s;
     536           0 :     y = shifti(y,-r);
     537             :   }
     538       34954 :   return gc_long(av, krouodd(x, y, s));
     539             : }
     540             : 
     541             : long
     542    96276421 : kross(long x, long y)
     543             : {
     544             :   ulong yu;
     545    96276421 :   long s = 1;
     546             : 
     547    96276421 :   if (y <= 0)
     548             :   {
     549       68943 :     if (y == 0) return (labs(x)==1);
     550       68915 :     yu = (ulong)-y; if (x < 0) s = -1;
     551             :   }
     552             :   else
     553    96207478 :     yu = (ulong)y;
     554    96276393 :   if (!odd(yu))
     555             :   {
     556             :     long r;
     557    23414185 :     if (!odd(x)) return 0;
     558    16495645 :     r = vals(yu); yu >>= r;
     559    16495645 :     if (odd(r) && ome(x)) s = -s;
     560             :   }
     561    89357853 :   x %= (long)yu; if (x < 0) x += yu;
     562    89357853 :   return krouu_s((ulong)x, yu, s);
     563             : }
     564             : 
     565             : long
     566    96403748 : krouu(ulong x, ulong y)
     567             : {
     568             :   long r;
     569    96403748 :   if (odd(y)) return krouu_s(x, y, 1);
     570        5339 :   if (!odd(x)) return 0;
     571        5527 :   r = vals(y); y >>= r;
     572        5527 :   return krouu_s(x, y, (odd(r) && ome(x))? -1: 1);
     573             : }
     574             : 
     575             : /*********************************************************************/
     576             : /**                          HILBERT SYMBOL                         **/
     577             : /*********************************************************************/
     578             : /* x,y are t_INT or t_REAL */
     579             : static long
     580        7329 : mphilbertoo(GEN x, GEN y)
     581             : {
     582        7329 :   long sx = signe(x), sy = signe(y);
     583        7329 :   if (!sx || !sy) return 0;
     584        7329 :   return (sx < 0 && sy < 0)? -1: 1;
     585             : }
     586             : 
     587             : long
     588      135331 : hilbertii(GEN x, GEN y, GEN p)
     589             : {
     590             :   pari_sp av;
     591             :   long oddvx, oddvy, z;
     592             : 
     593      135331 :   if (!p) return mphilbertoo(x,y);
     594      128023 :   if (is_pm1(p) || signe(p) < 0) pari_err_PRIME("hilbertii",p);
     595      128023 :   if (!signe(x) || !signe(y)) return 0;
     596      128002 :   av = avma;
     597      128002 :   oddvx = odd(Z_pvalrem(x,p,&x));
     598      128002 :   oddvy = odd(Z_pvalrem(y,p,&y));
     599             :   /* x, y are p-units, compute hilbert(x * p^oddvx, y * p^oddvy, p) */
     600      128002 :   if (absequaliu(p, 2))
     601             :   {
     602       12103 :     z = (Mod4(x) == 3 && Mod4(y) == 3)? -1: 1;
     603       12103 :     if (oddvx && gome(y)) z = -z;
     604       12103 :     if (oddvy && gome(x)) z = -z;
     605             :   }
     606             :   else
     607             :   {
     608      115899 :     z = (oddvx && oddvy && mod4(p) == 3)? -1: 1;
     609      115899 :     if (oddvx && kronecker(y,p) < 0) z = -z;
     610      115899 :     if (oddvy && kronecker(x,p) < 0) z = -z;
     611             :   }
     612      128002 :   return gc_long(av, z);
     613             : }
     614             : 
     615             : static void
     616         196 : err_prec(void) { pari_err_PREC("hilbert"); }
     617             : static void
     618         161 : err_p(GEN p, GEN q) { pari_err_MODULUS("hilbert", p,q); }
     619             : static void
     620          56 : err_oo(GEN p) { pari_err_MODULUS("hilbert", p, strtoGENstr("oo")); }
     621             : 
     622             : /* x t_INTMOD, *pp = prime or NULL [ unset, set it to x.mod ].
     623             :  * Return lift(x) provided it's p-adic accuracy is large enough to decide
     624             :  * hilbert()'s value [ problem at p = 2 ] */
     625             : static GEN
     626         420 : lift_intmod(GEN x, GEN *pp)
     627             : {
     628         420 :   GEN p = *pp, N = gel(x,1);
     629         420 :   x = gel(x,2);
     630         420 :   if (!p)
     631             :   {
     632         266 :     *pp = p = N;
     633         266 :     switch(itos_or_0(p))
     634             :     {
     635         126 :       case 2:
     636         126 :       case 4: err_prec();
     637             :     }
     638         140 :     return x;
     639             :   }
     640         154 :   if (!signe(p)) err_oo(N);
     641         112 :   if (absequaliu(p,2))
     642          42 :   { if (vali(N) <= 2) err_prec(); }
     643             :   else
     644          70 :   { if (!dvdii(N,p)) err_p(N,p); }
     645          28 :   if (!signe(x)) err_prec();
     646          21 :   return x;
     647             : }
     648             : /* x t_PADIC, *pp = prime or NULL [ unset, set it to x.p ].
     649             :  * Return lift(x)*p^(v(x) mod 2) provided it's p-adic accuracy is large enough
     650             :  * to decide hilbert()'s value [ problem at p = 2 ]*/
     651             : static GEN
     652         210 : lift_padic(GEN x, GEN *pp)
     653             : {
     654         210 :   GEN p = *pp, q = gel(x,2), y = gel(x,4);
     655         210 :   if (!p) *pp = p = q;
     656         147 :   else if (!equalii(p,q)) err_p(p, q);
     657         105 :   if (absequaliu(p,2) && precp(x) <= 2) err_prec();
     658          70 :   if (!signe(y)) err_prec();
     659          70 :   return odd(valp(x))? mulii(p,y): y;
     660             : }
     661             : 
     662             : long
     663       60403 : hilbert(GEN x, GEN y, GEN p)
     664             : {
     665       60403 :   pari_sp av = avma;
     666       60403 :   long tx = typ(x), ty = typ(y);
     667             : 
     668       60403 :   if (p && typ(p) != t_INT) pari_err_TYPE("hilbert",p);
     669       60403 :   if (tx == t_REAL)
     670             :   {
     671          77 :     if (p && signe(p)) err_oo(p);
     672          63 :     switch (ty)
     673             :     {
     674           7 :       case t_INT:
     675           7 :       case t_REAL: return mphilbertoo(x,y);
     676           0 :       case t_FRAC: return mphilbertoo(x,gel(y,1));
     677          56 :       default: pari_err_TYPE2("hilbert",x,y);
     678             :     }
     679             :   }
     680       60326 :   if (ty == t_REAL)
     681             :   {
     682          14 :     if (p && signe(p)) err_oo(p);
     683          14 :     switch (tx)
     684             :     {
     685          14 :       case t_INT:
     686          14 :       case t_REAL: return mphilbertoo(x,y);
     687           0 :       case t_FRAC: return mphilbertoo(gel(x,1),y);
     688           0 :       default: pari_err_TYPE2("hilbert",x,y);
     689             :     }
     690             :   }
     691       60312 :   if (tx == t_INTMOD) { x = lift_intmod(x, &p); tx = t_INT; }
     692       60109 :   if (ty == t_INTMOD) { y = lift_intmod(y, &p); ty = t_INT; }
     693             : 
     694       60053 :   if (tx == t_PADIC) { x = lift_padic(x, &p); tx = t_INT; }
     695       59990 :   if (ty == t_PADIC) { y = lift_padic(y, &p); ty = t_INT; }
     696             : 
     697       59913 :   if (tx == t_FRAC) { tx = t_INT; x = p? mulii(gel(x,1),gel(x,2)): gel(x,1); }
     698       59913 :   if (ty == t_FRAC) { ty = t_INT; y = p? mulii(gel(y,1),gel(y,2)): gel(y,1); }
     699             : 
     700       59913 :   if (tx != t_INT || ty != t_INT) pari_err_TYPE2("hilbert",x,y);
     701       59913 :   if (p && !signe(p)) p = NULL;
     702       59913 :   return gc_long(av, hilbertii(x,y,p));
     703             : }
     704             : 
     705             : /*******************************************************************/
     706             : /*                       SQUARE ROOT MODULO p                      */
     707             : /*******************************************************************/
     708             : static void
     709     2080236 : checkp(ulong q, ulong p)
     710     2080236 : { if (!q) pari_err_PRIME("Fl_nonsquare",utoipos(p)); }
     711             : /* p = 1 (mod 4) prime, return the first quadratic nonresidue, a prime */
     712             : static ulong
     713    11208283 : nonsquare1_Fl(ulong p)
     714             : {
     715             :   forprime_t S;
     716             :   ulong q;
     717    11208283 :   if ((p & 7UL) != 1) return 2UL;
     718     4166732 :   q = p % 3; if (q == 2) return 3UL;
     719     1299494 :   checkp(q, p);
     720     1309201 :   q = p % 5; if (q == 2 || q == 3) return 5UL;
     721      481969 :   checkp(q, p);
     722      481956 :   q = p % 7; if (q != 4 && q >= 3) return 7UL;
     723      177424 :   checkp(q, p);
     724             :   /* log^2(2^64) < 1968 is enough under GRH (and p^(1/4)log(p) without it)*/
     725      177478 :   u_forprime_init(&S, 11, 1967);
     726      289395 :   while ((q = u_forprime_next(&S)))
     727             :   {
     728      289392 :     if (krouu(q, p) < 0) return q;
     729      111911 :     checkp(q, p);
     730             :   }
     731           0 :   checkp(0, p);
     732             :   return 0; /*LCOV_EXCL_LINE*/
     733             : }
     734             : /* p > 2 a prime */
     735             : ulong
     736        7895 : nonsquare_Fl(ulong p)
     737        7895 : { return ((p & 3UL) == 3)? p-1: nonsquare1_Fl(p); }
     738             : 
     739             : /* allow pi = 0 */
     740             : ulong
     741      174704 : Fl_2gener_pre(ulong p, ulong pi)
     742             : {
     743      174704 :   ulong p1 = p-1;
     744      174704 :   long e = vals(p1);
     745      174680 :   if (e == 1) return p1;
     746       64904 :   return Fl_powu_pre(nonsquare1_Fl(p), p1 >> e, p, pi);
     747             : }
     748             : 
     749             : ulong
     750       67856 : Fl_2gener_pre_i(ulong  ns, ulong p, ulong pi)
     751             : {
     752       67856 :   ulong p1 = p-1;
     753       67856 :   long e = vals(p1);
     754       67856 :   if (e == 1) return p1;
     755       26020 :   return Fl_powu_pre(ns, p1 >> e, p, pi);
     756             : }
     757             : 
     758             : static ulong
     759    11856470 : Fl_sqrt_i(ulong a, ulong y, ulong p)
     760             : {
     761             :   long i, e, k;
     762             :   ulong p1, q, v, w;
     763             : 
     764    11856470 :   if (!a) return 0;
     765    10599180 :   p1 = p - 1; e = vals(p1);
     766    10599580 :   if (e == 0) /* p = 2 */
     767             :   {
     768      649926 :     if (p != 2) pari_err_PRIME("Fl_sqrt [modulus]",utoi(p));
     769      650976 :     return ((a & 1) == 0)? 0: 1;
     770             :   }
     771     9949654 :   if (e == 1)
     772             :   {
     773     4804480 :     v = Fl_powu(a, (p+1) >> 2, p);
     774     4804482 :     if (Fl_sqr(v, p) != a) return ~0UL;
     775     4799621 :     p1 = p - v; if (v > p1) v = p1;
     776     4799621 :     return v;
     777             :   }
     778     5145174 :   q = p1 >> e; /* q = (p-1)/2^oo is odd */
     779     5145174 :   p1 = Fl_powu(a, q >> 1, p); /* a ^ [(q-1)/2] */
     780     5145155 :   if (!p1) return 0;
     781     5145155 :   v = Fl_mul(a, p1, p);
     782     5145149 :   w = Fl_mul(v, p1, p);
     783     5145153 :   if (!y) y = Fl_powu(nonsquare1_Fl(p), q, p);
     784     8783191 :   while (w != 1)
     785             :   { /* a*w = v^2, y primitive 2^e-th root of 1
     786             :        a square --> w even power of y, hence w^(2^(e-1)) = 1 */
     787     3639871 :     p1 = Fl_sqr(w, p);
     788     5990641 :     for (k=1; p1 != 1 && k < e; k++) p1 = Fl_sqr(p1, p);
     789     3639884 :     if (k == e) return ~0UL;
     790             :     /* w ^ (2^k) = 1 --> w = y ^ (u * 2^(e-k)), u odd */
     791     3637807 :     p1 = y;
     792     4776728 :     for (i=1; i < e-k; i++) p1 = Fl_sqr(p1, p);
     793     3637807 :     y = Fl_sqr(p1, p); e = k;
     794     3637843 :     w = Fl_mul(y, w, p);
     795     3637853 :     v = Fl_mul(v, p1, p);
     796             :   }
     797     5143320 :   p1 = p - v; if (v > p1) v = p1;
     798     5143320 :   return v;
     799             : }
     800             : 
     801             : /* Tonelli-Shanks. Assume p is prime and (a,p) != -1. Allow pi = 0 */
     802             : ulong
     803    32979628 : Fl_sqrt_pre_i(ulong a, ulong y, ulong p, ulong pi)
     804             : {
     805             :   long i, e, k;
     806             :   ulong p1, q, v, w;
     807             : 
     808    32979628 :   if (!pi) return Fl_sqrt_i(a, y, p);
     809    21123349 :   if (!a) return 0;
     810    21000478 :   p1 = p - 1; e = vals(p1);
     811    21001224 :   if (e == 0) /* p = 2 */
     812             :   {
     813           0 :     if (p != 2) pari_err_PRIME("Fl_sqrt [modulus]",utoi(p));
     814           0 :     return ((a & 1) == 0)? 0: 1;
     815             :   }
     816    21015182 :   if (e == 1)
     817             :   {
     818    14958900 :     v = Fl_powu_pre(a, (p+1) >> 2, p, pi);
     819    14927711 :     if (Fl_sqr_pre(v, p, pi) != a) return ~0UL;
     820    14934854 :     p1 = p - v; if (v > p1) v = p1;
     821    14934854 :     return v;
     822             :   }
     823     6056282 :   q = p1 >> e; /* q = (p-1)/2^oo is odd */
     824     6056282 :   p1 = Fl_powu_pre(a, q >> 1, p, pi); /* a ^ [(q-1)/2] */
     825     6050003 :   if (!p1) return 0;
     826     6050003 :   v = Fl_mul_pre(a, p1, p, pi);
     827     6050310 :   w = Fl_mul_pre(v, p1, p, pi);
     828     6047468 :   if (!y) y = Fl_powu_pre(nonsquare1_Fl(p), q, p, pi);
     829    11527099 :   while (w != 1)
     830             :   { /* a*w = v^2, y primitive 2^e-th root of 1
     831             :        a square --> w even power of y, hence w^(2^(e-1)) = 1 */
     832     5476987 :     p1 = Fl_sqr_pre(w,p,pi);
     833    10245870 :     for (k=1; p1 != 1 && k < e; k++) p1 = Fl_sqr_pre(p1,p,pi);
     834     5477300 :     if (k == e) return ~0UL;
     835             :     /* w ^ (2^k) = 1 --> w = y ^ (u * 2^(e-k)), u odd */
     836     5477208 :     p1 = y;
     837     7180318 :     for (i=1; i < e-k; i++) p1 = Fl_sqr_pre(p1, p, pi);
     838     5477268 :     y = Fl_sqr_pre(p1, p, pi); e = k;
     839     5479018 :     w = Fl_mul_pre(y, w, p, pi);
     840     5477211 :     v = Fl_mul_pre(v, p1, p, pi);
     841             :   }
     842     6050112 :   p1 = p - v; if (v > p1) v = p1;
     843     6050112 :   return v;
     844             : }
     845             : 
     846             : ulong
     847    11871108 : Fl_sqrt(ulong a, ulong p)
     848    11871108 : { ulong pi = (p & HIGHMASK)? get_Fl_red(p): 0; return Fl_sqrt_pre_i(a, 0, p, pi); }
     849             : 
     850             : ulong
     851    20948129 : Fl_sqrt_pre(ulong a, ulong p, ulong pi)
     852    20948129 : { return Fl_sqrt_pre_i(a, 0, p, pi); }
     853             : 
     854             : /* allow pi = 0 */
     855             : static ulong
     856      140504 : Fl_lgener_pre_all(ulong l, long e, ulong r, ulong p, ulong pi, ulong *pt_m)
     857             : {
     858      140504 :   ulong x, y, m, le1 = upowuu(l, e-1);
     859      140504 :   for (x = 2; ; x++)
     860             :   {
     861      171263 :     y = Fl_powu_pre(x, r, p, pi);
     862      171263 :     if (y==1) continue;
     863      154237 :     m = Fl_powu_pre(y, le1, p, pi);
     864      154237 :     if (m != 1) break;
     865             :   }
     866      140504 :   *pt_m = m; return y;
     867             : }
     868             : 
     869             : /* solve x^l = a , l prime in G of order q.
     870             :  *
     871             :  * q =  (l^e)*r, e >= 1, (r,l) = 1
     872             :  * y generates the l-Sylow of G
     873             :  * m = y^(l^(e-1)) != 1 */
     874             : static ulong
     875      223736 : Fl_sqrtl_raw(ulong a, ulong l, ulong e, ulong r, ulong p, ulong pi, ulong y, ulong m)
     876             : {
     877             :   ulong u2, p1, v, w, z, dl;
     878      223736 :   if (a==0) return a;
     879      223730 :   u2 = Fl_inv(l%r, r);
     880      223732 :   v = Fl_powu_pre(a, u2, p, pi);
     881      223734 :   w = Fl_powu_pre(v, l, p, pi);
     882      223735 :   w = pi? Fl_mul_pre(w, Fl_inv(a, p), p, pi): Fl_div(w, a, p);
     883      223719 :   if (w==1) return v;
     884      138423 :   if (y==0) y = Fl_lgener_pre_all(l, e, r, p, pi, &m);
     885      163380 :   while (w!=1)
     886             :   {
     887      143767 :     ulong k = 0;
     888      143767 :     p1 = w;
     889             :     do
     890             :     {
     891      187017 :       z = p1; p1 = Fl_powu_pre(p1, l, p, pi);
     892      187018 :       if (++k == e) return ULONG_MAX;
     893       68207 :     } while (p1!=1);
     894       24957 :     dl = Fl_log_pre(z, m, l, p, pi);
     895       24957 :     dl = Fl_neg(dl, l);
     896       24957 :     p1 = Fl_powu_pre(y,dl*upowuu(l,e-k-1),p,pi);
     897       24956 :     m = Fl_powu_pre(m, dl, p, pi);
     898       24956 :     e = k;
     899       24956 :     v = pi? Fl_mul_pre(p1,v,p,pi): Fl_mul(p1,v,p);
     900       24956 :     y = Fl_powu_pre(p1,l,p,pi);
     901       24957 :     w = pi? Fl_mul_pre(y,w,p,pi): Fl_mul(y,w,p);
     902             :   }
     903       19613 :   return v;
     904             : }
     905             : 
     906             : /* allow pi = 0 */
     907             : static ulong
     908      221369 : Fl_sqrtl_i(ulong a, ulong l, ulong p, ulong pi, ulong y, ulong m)
     909             : {
     910      221369 :   ulong r, e = u_lvalrem(p-1, l, &r);
     911      221369 :   return Fl_sqrtl_raw(a, l, e, r, p, pi, y, m);
     912             : }
     913             : /* allow pi = 0 */
     914             : ulong
     915      221369 : Fl_sqrtl_pre(ulong a, ulong l, ulong p, ulong pi)
     916      221369 : { return Fl_sqrtl_i(a, l, p, pi, 0, 0); }
     917             : 
     918             : ulong
     919           0 : Fl_sqrtl(ulong a, ulong l, ulong p)
     920           0 : { ulong pi = (p & HIGHMASK)? get_Fl_red(p): 0;
     921           0 :   return Fl_sqrtl_i(a, l, p, pi, 0, 0); }
     922             : 
     923             : /* allow pi = 0 */
     924             : ulong
     925      196266 : Fl_sqrtn_pre(ulong a, long n, ulong p, ulong pi, ulong *zetan)
     926             : {
     927      196266 :   ulong m, q = p-1, z;
     928      196266 :   ulong nn = n >= 0 ? (ulong)n: -(ulong)n;
     929      196266 :   if (a==0)
     930             :   {
     931      116389 :     if (n < 0) pari_err_INV("Fl_sqrtn", mkintmod(gen_0,utoi(p)));
     932      116382 :     if (zetan) *zetan = 1UL;
     933      116382 :     return 0;
     934             :   }
     935       79877 :   if (n==1)
     936             :   {
     937         420 :     if (zetan) *zetan = 1;
     938         420 :     return n < 0? Fl_inv(a,p): a;
     939             :   }
     940       79457 :   if (n==2)
     941             :   {
     942        6483 :     if (zetan) *zetan = p-1;
     943        6483 :     return Fl_sqrt_pre_i(a, 0, p, pi);
     944             :   }
     945       72974 :   if (a == 1 && !zetan) return a;
     946       43551 :   m = ugcd(nn, q);
     947       43551 :   z = 1;
     948       43551 :   if (m!=1)
     949             :   {
     950        2076 :     GEN F = factoru(m);
     951             :     long i, j, e;
     952             :     ulong r, zeta, y, l;
     953        4535 :     for (i = nbrows(F); i; i--)
     954             :     {
     955        2522 :       l = ucoeff(F,i,1);
     956        2522 :       j = ucoeff(F,i,2);
     957        2522 :       e = u_lvalrem(q,l, &r);
     958        2522 :       y = Fl_lgener_pre_all(l, e, r, p, pi, &zeta);
     959        2522 :       if (zetan)
     960             :       {
     961         770 :         ulong Y = Fl_powu_pre(y, upowuu(l,e-j), p, pi);
     962         770 :         z = pi? Fl_mul_pre(z, Y, p, pi): Fl_mul(z, Y, p);
     963             :       }
     964        2522 :       if (a!=1)
     965             :         do
     966             :         {
     967        2367 :           a = Fl_sqrtl_raw(a, l, e, r, p, pi, y, zeta);
     968        2353 :           if (a==ULONG_MAX) return ULONG_MAX;
     969        2304 :         } while (--j);
     970             :     }
     971             :   }
     972       43488 :   if (m != nn)
     973             :   {
     974       41496 :     ulong qm = q/m, nm = (nn/m) % qm;
     975       41496 :     a = Fl_powu_pre(a, Fl_inv(nm, qm), p, pi);
     976             :   }
     977       43488 :   if (n < 0) a = Fl_inv(a, p);
     978       43488 :   if (zetan) *zetan = z;
     979       43488 :   return a;
     980             : }
     981             : 
     982             : ulong
     983      196266 : Fl_sqrtn(ulong a, long n, ulong p, ulong *zetan)
     984             : {
     985      196266 :   ulong pi = (p & HIGHMASK)? get_Fl_red(p): 0;
     986      196266 :   return Fl_sqrtn_pre(a, n, p, pi, zetan);
     987             : }
     988             : 
     989             : /* Cipolla is better than Tonelli-Shanks when e = v_2(p-1) is "too big".
     990             :  * Otherwise, is a constant times worse; for p = 3 (mod 4), is about 3 times worse,
     991             :  * and in average is about 2 or 2.5 times worse. But try both algorithms for
     992             :  * S(n) = (2^n+3)^2-8 with n = 750, 771, 779, 790, 874, 1176, 1728, 2604, etc.
     993             :  *
     994             :  * If X^2 := t^2 - a  is not a square in F_p (so X is in F_p^2), then
     995             :  *   (t+X)^(p+1) = (t-X)(t+X) = a,   hence  sqrt(a) = (t+X)^((p+1)/2)  in F_p^2.
     996             :  * If (a|p)=1, then sqrt(a) is in F_p.
     997             :  * cf: LNCS 2286, pp 430-434 (2002)  [Gonzalo Tornaria] */
     998             : 
     999             : /* compute y^2, y = y[1] + y[2] X */
    1000             : static GEN
    1001           0 : sqrt_Cipolla_sqr(void *data, GEN y)
    1002             : {
    1003           0 :   GEN u = gel(y,1), v = gel(y,2), p = gel(data,2), n = gel(data,3);
    1004           0 :   GEN u2 = sqri(u), v2 = sqri(v);
    1005           0 :   v = subii(sqri(addii(v,u)), addii(u2,v2));
    1006           0 :   u = addii(u2, mulii(v2,n));
    1007           0 :   retmkvec2(modii(u,p), modii(v,p));
    1008             : }
    1009             : /* compute (t+X) y^2 */
    1010             : static GEN
    1011           0 : sqrt_Cipolla_msqr(void *data, GEN y)
    1012             : {
    1013           0 :   GEN u = gel(y,1), v = gel(y,2), a = gel(data,1), p = gel(data,2);
    1014           0 :   ulong t = gel(data,4)[2];
    1015           0 :   GEN d = addii(u, mului(t,v)), d2 = sqri(d);
    1016           0 :   GEN b = remii(mulii(a,v), p);
    1017           0 :   u = subii(mului(t,d2), mulii(b,addii(u,d)));
    1018           0 :   v = subii(d2, mulii(b,v));
    1019           0 :   retmkvec2(modii(u,p), modii(v,p));
    1020             : }
    1021             : /* assume a reduced mod p [ otherwise correct but inefficient ] */
    1022             : static GEN
    1023           0 : sqrt_Cipolla(GEN a, GEN p)
    1024             : {
    1025             :   pari_sp av;
    1026             :   GEN u, n, y, pov2;
    1027             :   ulong t;
    1028             : 
    1029           0 :   if (kronecker(a, p) < 0) return NULL;
    1030           0 :   pov2 = shifti(p,-1); /* center to avoid multiplying by huge base*/
    1031           0 :   if (cmpii(a,pov2) > 0) a = subii(a,p);
    1032           0 :   av = avma;
    1033           0 :   for (t=1; ; t++, set_avma(av))
    1034             :   {
    1035           0 :     n = subsi((long)(t*t), a);
    1036           0 :     if (kronecker(n, p) < 0) break;
    1037             :   }
    1038             : 
    1039             :   /* compute (t+X)^((p-1)/2) =: u+vX */
    1040           0 :   u = utoipos(t);
    1041           0 :   y = gen_pow_fold(mkvec2(u, gen_1), pov2, mkvec4(a,p,n,u),
    1042             :                    sqrt_Cipolla_sqr, sqrt_Cipolla_msqr);
    1043             :   /* Now u+vX = (t+X)^((p-1)/2); thus
    1044             :    *   (u+vX)(t+X) = sqrt(a) + 0 X
    1045             :    * Whence,
    1046             :    *   sqrt(a) = (u+vt)t - v*a
    1047             :    *   0       = (u+vt)
    1048             :    * Thus a square root is v*a */
    1049           0 :   return Fp_mul(gel(y,2), a, p);
    1050             : }
    1051             : 
    1052             : /* Return NULL if p is found to be composite.
    1053             :  * p odd, q = (p-1)/2^oo is odd */
    1054             : static GEN
    1055        5904 : Fp_2gener_all(GEN q, GEN p)
    1056             : {
    1057             :   long k;
    1058        5904 :   for (k = 2;; k++)
    1059       11607 :   {
    1060       17511 :     long i = kroui(k, p);
    1061       17511 :     if (i < 0) return Fp_pow(utoipos(k), q, p);
    1062       11607 :     if (i == 0) return NULL;
    1063             :   }
    1064             : }
    1065             : 
    1066             : /* Return NULL if p is found to be composite */
    1067             : GEN
    1068        3192 : Fp_2gener(GEN p)
    1069             : {
    1070        3192 :   GEN q = subiu(p, 1);
    1071        3192 :   long e = Z_lvalrem(q, 2, &q);
    1072        3192 :   if (e == 0 && !equaliu(p,2)) return NULL;
    1073        3192 :   return Fp_2gener_all(q, p);
    1074             : }
    1075             : 
    1076             : GEN
    1077       19931 : Fp_2gener_i(GEN ns, GEN p)
    1078             : {
    1079       19931 :   GEN q = subiu(p,1);
    1080       19931 :   long e = vali(q);
    1081       19931 :   if (e == 1) return q;
    1082       18546 :   return Fp_pow(ns, shifti(q,-e), p);
    1083             : }
    1084             : 
    1085             : static GEN
    1086        1472 : nonsquare_Fp(GEN p)
    1087             : {
    1088             :   forprime_t T;
    1089             :   ulong a;
    1090        1472 :   if (mod4(p)==3) return gen_m1;
    1091        1472 :   if (mod8(p)==5) return gen_2;
    1092         721 :   u_forprime_init(&T, 3, ULONG_MAX);
    1093        1382 :   while((a = u_forprime_next(&T)))
    1094        1382 :     if (kroui(a,p) < 0) return utoi(a);
    1095           0 :   pari_err_PRIME("Fp_sqrt [modulus]",p);
    1096             :   return NULL; /* LCOV_EXCL_LINE */
    1097             : }
    1098             : 
    1099             : static GEN
    1100         796 : Fp_rootsof1(ulong l, GEN p)
    1101             : {
    1102         796 :   GEN z, pl = diviuexact(subis(p,1),l);
    1103             :   ulong a;
    1104             :   forprime_t T;
    1105         796 :   u_forprime_init(&T, 3, ULONG_MAX);
    1106        1062 :   while((a = u_forprime_next(&T)))
    1107             :   {
    1108        1062 :     z = Fp_pow(utoi(a), pl, p);
    1109        1062 :     if (!equali1(z)) return z;
    1110             :   }
    1111           0 :   pari_err_PRIME("Fp_sqrt [modulus]",p);
    1112             :   return NULL; /* LCOV_EXCL_LINE */
    1113             : }
    1114             : 
    1115             : static GEN
    1116         334 : Fp_gausssum(long D, GEN p)
    1117             : {
    1118         334 :   long i, l = labs(D);
    1119         334 :   GEN z = Fp_rootsof1(l, p);
    1120         334 :   GEN s = z, x = z;
    1121        3020 :   for(i = 2; i < l; i++)
    1122             :   {
    1123        2686 :     long k = kross(i,l);
    1124        2686 :     x = mulii(x, z);
    1125        2686 :     if (k==1) s = addii(s, x);
    1126        1510 :     else if (k==-1) s = subii(s, x);
    1127             :   }
    1128         334 :   return s;
    1129             : }
    1130             : 
    1131             : static GEN
    1132       19199 : Fp_sqrts(long a, GEN p)
    1133             : {
    1134       19199 :   long v = vals(a)>>1;
    1135       19199 :   GEN r = gen_0;
    1136       19199 :   a >>= v << 1;
    1137       19199 :   switch(a)
    1138             :   {
    1139           8 :     case 1:
    1140           8 :       r = gen_1;
    1141           8 :       break;
    1142        1128 :     case -1:
    1143        1128 :       if (mod4(p)==1)
    1144        1128 :         r = Fp_pow(nonsquare_Fp(p), shifti(p,-2),p);
    1145             :       else
    1146           0 :         r = NULL;
    1147        1128 :       break;
    1148         140 :     case 2:
    1149         140 :       if (mod8(p)==1)
    1150             :       {
    1151         140 :         GEN z = Fp_pow(nonsquare_Fp(p), shifti(p,-3),p);
    1152         140 :         r = Fp_mul(z,Fp_sub(gen_1,Fp_sqr(z,p),p),p);
    1153           0 :       } else if (mod8(p)==7)
    1154           0 :         r = Fp_pow(gen_2, shifti(addiu(p,1),-2),p);
    1155             :       else
    1156           0 :         return NULL;
    1157         140 :       break;
    1158         204 :     case -2:
    1159         204 :       if (mod8(p)==1)
    1160             :       {
    1161         204 :         GEN z = Fp_pow(nonsquare_Fp(p), shifti(p,-3),p);
    1162         204 :         r = Fp_mul(z,Fp_add(gen_1,Fp_sqr(z,p),p),p);
    1163           0 :       } else if (mod8(p)==3)
    1164           0 :         r = Fp_pow(gen_m2, shifti(addiu(p,1),-2),p);
    1165             :       else
    1166           0 :         return NULL;
    1167         204 :       break;
    1168         462 :     case -3:
    1169         462 :       if (umodiu(p,3)==1)
    1170             :       {
    1171         462 :         GEN z = Fp_rootsof1(3, p);
    1172         462 :         r = Fp_sub(z,Fp_sqr(z,p),p);
    1173             :       }
    1174             :       else
    1175           0 :         return NULL;
    1176         462 :       break;
    1177        2213 :     case 5: case 13: case 17: case 21: case 29: case 33:
    1178             :     case -7: case -11: case -15: case -19: case -23:
    1179        2213 :       if (umodiu(p,labs(a))==1)
    1180         334 :         r = Fp_gausssum(a,p);
    1181             :       else
    1182        1880 :         return gen_0;
    1183         334 :       break;
    1184       15044 :     default:
    1185       15044 :       return gen_0;
    1186             :   }
    1187        2276 :   return remii(shifti(r, v), p);
    1188             : }
    1189             : 
    1190             : static GEN
    1191       77468 : Fp_sqrt_ii(GEN a, GEN y, GEN p)
    1192             : {
    1193       77468 :   pari_sp av = avma;
    1194       77468 :   GEN  q, v, w, p1 = subiu(p,1);
    1195       77472 :   long i, k, e = vali(p1), as;
    1196             : 
    1197             :   /* direct formulas more efficient */
    1198       77479 :   if (e == 0) pari_err_PRIME("Fp_sqrt [modulus]",p); /* p != 2 */
    1199       77480 :   if (e == 1)
    1200             :   {
    1201       19091 :     q = addiu(shifti(p1,-2),1); /* (p+1) / 4 */
    1202       19091 :     v = Fp_pow(a, q, p);
    1203             :     /* must check equality in case (a/p) = -1 or p not prime */
    1204       19099 :     av = avma; e = equalii(Fp_sqr(v,p), a); set_avma(av);
    1205       19098 :     return e? v: NULL;
    1206             :   }
    1207       58389 :   as = itos_or_0(a);
    1208       58389 :   if (!as) as = itos_or_0(subii(a,p));
    1209       58394 :   if (as)
    1210             :   {
    1211       19199 :     GEN res = Fp_sqrts(as, p);
    1212       19199 :     if (!res) return gc_NULL(av);
    1213       19199 :     if (signe(res)) return gerepileupto(av, res);
    1214             :   }
    1215       56118 :   if (e == 2)
    1216             :   { /* Atkin's formula */
    1217       17286 :     GEN I, a2 = shifti(a,1);
    1218       17284 :     if (cmpii(a2,p) >= 0) a2 = subii(a2,p);
    1219       17284 :     q = shifti(p1, -3); /* (p-5)/8 */
    1220       17282 :     v = Fp_pow(a2, q, p);
    1221       17288 :     I = Fp_mul(a2, Fp_sqr(v,p), p); /* I^2 = -1 */
    1222       17288 :     v = Fp_mul(a, Fp_mul(v, subiu(I,1), p), p);
    1223             :     /* must check equality in case (a/p) = -1 or p not prime */
    1224       17288 :     av = avma; e = equalii(Fp_sqr(v,p), a); set_avma(av);
    1225       17288 :     return e? v: NULL;
    1226             :   }
    1227             :   /* On average, Cipolla is better than Tonelli/Shanks if and only if
    1228             :    * e(e-1) > 8*log2(n)+20, see LNCS 2286 pp 430 [GTL] */
    1229       38832 :   if (e*(e-1) > 20 + 8 * expi(p)) return sqrt_Cipolla(a,p);
    1230             :   /* Tonelli-Shanks */
    1231       38832 :   av = avma; q = shifti(p1,-e); /* q = (p-1)/2^oo is odd */
    1232       38832 :   if (!y)
    1233             :   {
    1234        2712 :     y = Fp_2gener_all(q, p);
    1235        2712 :     if (!y) pari_err_PRIME("Fp_sqrt [modulus]",p);
    1236             :   }
    1237       38832 :   p1 = Fp_pow(a, shifti(q,-1), p); /* a ^ (q-1)/2 */
    1238       38832 :   v = Fp_mul(a, p1, p);
    1239       38832 :   w = Fp_mul(v, p1, p);
    1240       92603 :   while (!equali1(w))
    1241             :   { /* a*w = v^2, y primitive 2^e-th root of 1
    1242             :        a square --> w even power of y, hence w^(2^(e-1)) = 1 */
    1243       53814 :     p1 = Fp_sqr(w,p);
    1244      110469 :     for (k=1; !equali1(p1) && k < e; k++) p1 = Fp_sqr(p1,p);
    1245       53811 :     if (k == e) return NULL; /* p composite or (a/p) != 1 */
    1246             :     /* w ^ (2^k) = 1 --> w = y ^ (u * 2^(e-k)), u odd */
    1247       53768 :     p1 = y;
    1248       76499 :     for (i=1; i < e-k; i++) p1 = Fp_sqr(p1,p);
    1249       53768 :     y = Fp_sqr(p1, p); e = k;
    1250       53773 :     w = Fp_mul(y, w, p);
    1251       53772 :     v = Fp_mul(v, p1, p);
    1252       53772 :     if (gc_needed(av,1))
    1253             :     {
    1254           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"Fp_sqrt");
    1255           0 :       gerepileall(av,3, &y,&w,&v);
    1256             :     }
    1257             :   }
    1258       38789 :   return v;
    1259             : }
    1260             : 
    1261             : /* Assume p is prime and return NULL if (a,p) = -1; y = NULL or generator
    1262             :  * of Fp^* 2-Sylow */
    1263             : GEN
    1264     4597904 : Fp_sqrt_i(GEN a, GEN y, GEN p)
    1265             : {
    1266     4597904 :   pari_sp av = avma, av2;
    1267             :   GEN q;
    1268             : 
    1269     4597904 :   if (lgefint(p) == 3)
    1270             :   {
    1271     4520321 :     ulong pp = uel(p,2), u = umodiu(a, pp);
    1272     4520337 :     if (!u) return gen_0;
    1273     3309176 :     u = Fl_sqrt(u, pp);
    1274     3309306 :     return (u == ~0UL)? NULL: utoipos(u);
    1275             :   }
    1276       77583 :   a = modii(a, p); if (!signe(a)) return gen_0;
    1277       77468 :   a = Fp_sqrt_ii(a, y, p); if (!a) return gc_NULL(av);
    1278             :   /* smallest square root */
    1279       77075 :   av2 = avma; q = subii(p, a);
    1280       77075 :   if (cmpii(a, q) > 0) a = q; else set_avma(av2);
    1281       77076 :   return gerepileuptoint(av, a);
    1282             : }
    1283             : GEN
    1284     4540894 : Fp_sqrt(GEN a, GEN p) { return Fp_sqrt_i(a, NULL, p); }
    1285             : 
    1286             : /*********************************************************************/
    1287             : /**                        GCD & BEZOUT                             **/
    1288             : /*********************************************************************/
    1289             : 
    1290             : GEN
    1291    54042181 : lcmii(GEN x, GEN y)
    1292             : {
    1293             :   pari_sp av;
    1294             :   GEN a, b;
    1295    54042181 :   if (!signe(x) || !signe(y)) return gen_0;
    1296    54042185 :   av = avma; a = gcdii(x,y);
    1297    54041027 :   if (absequalii(a,y)) { set_avma(av); return absi(x); }
    1298    12123541 :   if (!equali1(a)) y = diviiexact(y,a);
    1299    12123532 :   b = mulii(x,y); setabssign(b); return gerepileuptoint(av, b);
    1300             : }
    1301             : 
    1302             : /* given x in assume 0 < x < N; return u in (Z/NZ)^* such that u x = gcd(x,N) (mod N);
    1303             :  * set *pd = gcd(x,N) */
    1304             : GEN
    1305     5803968 : Fp_invgen(GEN x, GEN N, GEN *pd)
    1306             : {
    1307             :   GEN d, d0, e, v;
    1308     5803968 :   if (lgefint(N) == 3)
    1309             :   {
    1310     5047387 :     ulong dd, NN = N[2], xx = umodiu(x,NN);
    1311     5047544 :     if (!xx) { *pd = N; return gen_0; }
    1312     5047544 :     xx = Fl_invgen(xx, NN, &dd);
    1313     5048812 :     *pd = utoi(dd); return utoi(xx);
    1314             :   }
    1315      756581 :   *pd = d = bezout(x, N, &v, NULL);
    1316      756590 :   if (equali1(d)) return v;
    1317             :   /* vx = gcd(x,N) (mod N), v coprime to N/d but need not be coprime to N */
    1318      662092 :   e = diviiexact(N,d);
    1319      662092 :   d0 = Z_ppo(d, e); /* d = d0 d1, d0 coprime to N/d, rad(d1) | N/d */
    1320      662092 :   if (equali1(d0)) return v;
    1321      524809 :   if (!equalii(d,d0)) e = lcmii(e, diviiexact(d,d0));
    1322      524809 :   return Z_chinese_coprime(v, gen_1, e, d0, mulii(e,d0));
    1323             : }
    1324             : 
    1325             : /*********************************************************************/
    1326             : /**                      CHINESE REMAINDERS                         **/
    1327             : /*********************************************************************/
    1328             : 
    1329             : /* Chinese Remainder Theorem.  x and y must have the same type (integermod,
    1330             :  * polymod, or polynomial/vector/matrix recursively constructed with these
    1331             :  * as coefficients). Creates (with the same type) a z in the same residue
    1332             :  * class as x and the same residue class as y, if it is possible.
    1333             :  *
    1334             :  * We also allow (during recursion) two identical objects even if they are
    1335             :  * not integermod or polymod. For example:
    1336             :  *
    1337             :  * ? x = [1, Mod(5, 11), Mod(X + Mod(2, 7), X^2 + 1)];
    1338             :  * ? y = [1, Mod(7, 17), Mod(X + Mod(0, 3), X^2 + 1)];
    1339             :  * ? chinese(x, y)
    1340             :  * %3 = [1, Mod(16, 187), Mod(X + mod(9, 21), X^2 + 1)] */
    1341             : 
    1342             : static GEN
    1343     2410183 : gen_chinese(GEN x, GEN(*f)(GEN,GEN))
    1344             : {
    1345     2410183 :   GEN z = gassoc_proto(f,x,NULL);
    1346     2410176 :   if (z == gen_1) retmkintmod(gen_0,gen_1);
    1347     2410141 :   return z;
    1348             : }
    1349             : 
    1350             : /* x t_INTMOD, y t_POLMOD; promote x to t_POLMOD mod Pol(x.mod) then
    1351             :  * call chinese: makes Mod(0,1) a better "neutral" element */
    1352             : static GEN
    1353          21 : chinese_intpol(GEN x,GEN y)
    1354             : {
    1355          21 :   pari_sp av = avma;
    1356          21 :   GEN z = mkpolmod(gel(x,2), scalarpol_shallow(gel(x,1), varn(gel(y,1))));
    1357          21 :   return gerepileupto(av, chinese(z, y));
    1358             : }
    1359             : 
    1360             : GEN
    1361        2345 : chinese1(GEN x) { return gen_chinese(x,chinese); }
    1362             : 
    1363             : GEN
    1364        5425 : chinese(GEN x, GEN y)
    1365             : {
    1366             :   pari_sp av;
    1367        5425 :   long tx = typ(x), ty;
    1368             :   GEN z,p1,p2,d,u,v;
    1369             : 
    1370        5425 :   if (!y) return chinese1(x);
    1371        5376 :   if (gidentical(x,y)) return gcopy(x);
    1372        5369 :   ty = typ(y);
    1373        5369 :   if (tx == ty) switch(tx)
    1374             :   {
    1375        3822 :     case t_POLMOD:
    1376             :     {
    1377        3822 :       GEN A = gel(x,1), B = gel(y,1);
    1378        3822 :       GEN a = gel(x,2), b = gel(y,2);
    1379        3822 :       if (varn(A)!=varn(B)) pari_err_VAR("chinese",A,B);
    1380        3822 :       if (RgX_equal(A,B)) retmkpolmod(chinese(a,b), gcopy(A)); /*same modulus*/
    1381        3822 :       av = avma;
    1382        3822 :       d = RgX_extgcd(A,B,&u,&v);
    1383        3822 :       p2 = gsub(b, a);
    1384        3822 :       if (!gequal0(gmod(p2, d))) break;
    1385        3822 :       p1 = gdiv(A,d);
    1386        3822 :       p2 = gadd(a, gmul(gmul(u,p1), p2));
    1387             : 
    1388        3822 :       z = cgetg(3, t_POLMOD);
    1389        3822 :       gel(z,1) = gmul(p1,B);
    1390        3822 :       gel(z,2) = gmod(p2,gel(z,1));
    1391        3822 :       return gerepileupto(av, z);
    1392             :     }
    1393        1505 :     case t_INTMOD:
    1394             :     {
    1395        1505 :       GEN A = gel(x,1), B = gel(y,1);
    1396        1505 :       GEN a = gel(x,2), b = gel(y,2), c, d, C, U;
    1397        1505 :       z = cgetg(3,t_INTMOD);
    1398        1505 :       Z_chinese_pre(A, B, &C, &U, &d);
    1399        1505 :       c = Z_chinese_post(a, b, C, U, d);
    1400        1505 :       if (!c) pari_err_OP("chinese", x,y);
    1401        1505 :       set_avma((pari_sp)z);
    1402        1505 :       gel(z,1) = icopy(C);
    1403        1505 :       gel(z,2) = icopy(c); return z;
    1404             :     }
    1405          14 :     case t_POL:
    1406             :     {
    1407          14 :       long i, lx = lg(x), ly = lg(y);
    1408          14 :       if (varn(x) != varn(y)) break;
    1409          14 :       if (lx < ly) { swap(x,y); lswap(lx,ly); }
    1410          14 :       z = cgetg(lx, t_POL); z[1] = x[1];
    1411          42 :       for (i=2; i<ly; i++) gel(z,i) = chinese(gel(x,i),gel(y,i));
    1412          14 :       if (i < lx)
    1413             :       {
    1414          14 :         GEN _0 = Rg_get_0(y);
    1415          28 :         for (   ; i<lx; i++) gel(z,i) = chinese(gel(x,i),_0);
    1416             :       }
    1417          14 :       return z;
    1418             :     }
    1419           7 :     case t_VEC: case t_COL: case t_MAT:
    1420             :     {
    1421             :       long i, lx;
    1422           7 :       z = cgetg_copy(x, &lx); if (lx!=lg(y)) break;
    1423          21 :       for (i=1; i<lx; i++) gel(z,i) = chinese(gel(x,i),gel(y,i));
    1424           7 :       return z;
    1425             :     }
    1426             :   }
    1427          21 :   if (tx == t_POLMOD && ty == t_INTMOD) return chinese_intpol(y,x);
    1428           7 :   if (ty == t_POLMOD && tx == t_INTMOD) return chinese_intpol(x,y);
    1429           0 :   pari_err_OP("chinese",x,y);
    1430             :   return NULL; /* LCOV_EXCL_LINE */
    1431             : }
    1432             : 
    1433             : /* init chinese(Mod(.,A), Mod(.,B)) */
    1434             : void
    1435      271318 : Z_chinese_pre(GEN A, GEN B, GEN *pC, GEN *pU, GEN *pd)
    1436             : {
    1437      271318 :   GEN u, d = bezout(A,B,&u,NULL); /* U = u(A/d), u(A/d) + v(B/d) = 1 */
    1438      271319 :   GEN t = diviiexact(A,d);
    1439      271306 :   *pU = mulii(u, t);
    1440      271303 :   *pC = mulii(t, B);
    1441      271308 :   if (pd) *pd = d;
    1442      271308 : }
    1443             : /* Assume C = lcm(A, B), U = 0 mod (A/d), U = 1 mod (B/d), a = b mod d,
    1444             :  * where d = gcd(A,B) or NULL, return x = a (mod A), b (mod B).
    1445             :  * If d not NULL, check whether a = b mod d. */
    1446             : GEN
    1447     2986204 : Z_chinese_post(GEN a, GEN b, GEN C, GEN U, GEN d)
    1448             : {
    1449             :   GEN b_a;
    1450     2986204 :   if (!signe(a))
    1451             :   {
    1452      794810 :     if (d && !dvdii(b, d)) return NULL;
    1453      794810 :     return Fp_mul(b, U, C);
    1454             :   }
    1455     2191394 :   b_a = subii(b,a);
    1456     2191394 :   if (d && !dvdii(b_a, d)) return NULL;
    1457     2191394 :   return modii(addii(a, mulii(U, b_a)), C);
    1458             : }
    1459             : static ulong
    1460     1523739 : u_chinese_post(ulong a, ulong b, ulong C, ulong U)
    1461             : {
    1462     1523739 :   if (!a) return Fl_mul(b, U, C);
    1463     1523739 :   return Fl_add(a, Fl_mul(U, Fl_sub(b,a,C), C), C);
    1464             : }
    1465             : 
    1466             : GEN
    1467        2142 : Z_chinese(GEN a, GEN b, GEN A, GEN B)
    1468             : {
    1469        2142 :   pari_sp av = avma;
    1470        2142 :   GEN C, U; Z_chinese_pre(A, B, &C, &U, NULL);
    1471        2142 :   return gerepileuptoint(av, Z_chinese_post(a,b, C, U, NULL));
    1472             : }
    1473             : GEN
    1474      267615 : Z_chinese_all(GEN a, GEN b, GEN A, GEN B, GEN *pC)
    1475             : {
    1476      267615 :   GEN U; Z_chinese_pre(A, B, pC, &U, NULL);
    1477      267605 :   return Z_chinese_post(a,b, *pC, U, NULL);
    1478             : }
    1479             : 
    1480             : /* return lift(chinese(a mod A, b mod B))
    1481             :  * assume(A,B)=1, a,b,A,B integers and C = A*B */
    1482             : GEN
    1483      526069 : Z_chinese_coprime(GEN a, GEN b, GEN A, GEN B, GEN C)
    1484             : {
    1485      526069 :   pari_sp av = avma;
    1486      526069 :   GEN U = mulii(Fp_inv(A,B), A);
    1487      526069 :   return gerepileuptoint(av, Z_chinese_post(a,b,C,U, NULL));
    1488             : }
    1489             : ulong
    1490     1523738 : u_chinese_coprime(ulong a, ulong b, ulong A, ulong B, ulong C)
    1491     1523738 : { return u_chinese_post(a,b,C, A * Fl_inv(A % B,B)); }
    1492             : 
    1493             : /* chinese1 for coprime moduli in Z */
    1494             : static GEN
    1495     2188561 : chinese1_coprime_Z_aux(GEN x, GEN y)
    1496             : {
    1497     2188561 :   GEN z = cgetg(3, t_INTMOD);
    1498     2188561 :   GEN A = gel(x,1), a = gel(x, 2);
    1499     2188561 :   GEN B = gel(y,1), b = gel(y, 2), C = mulii(A,B);
    1500     2188561 :   pari_sp av = avma;
    1501     2188561 :   GEN U = mulii(Fp_inv(A,B), A);
    1502     2188561 :   gel(z,2) = gerepileuptoint(av, Z_chinese_post(a,b,C,U, NULL));
    1503     2188561 :   gel(z,1) = C; return z;
    1504             : }
    1505             : GEN
    1506     2407838 : chinese1_coprime_Z(GEN x) {return gen_chinese(x,chinese1_coprime_Z_aux);}
    1507             : 
    1508             : /*********************************************************************/
    1509             : /**                    MODULAR EXPONENTIATION                       **/
    1510             : /*********************************************************************/
    1511             : /* xa ZV or nv */
    1512             : GEN
    1513     2428286 : ZV_producttree(GEN xa)
    1514             : {
    1515     2428286 :   long n = lg(xa)-1;
    1516     2428286 :   long m = n==1 ? 1: expu(n-1)+1;
    1517     2428286 :   GEN T = cgetg(m+1, t_VEC), t;
    1518             :   long i, j, k;
    1519     2428281 :   t = cgetg(((n+1)>>1)+1, t_VEC);
    1520     2428282 :   if (typ(xa)==t_VECSMALL)
    1521             :   {
    1522     3198379 :     for (j=1, k=1; k<n; j++, k+=2)
    1523     2052128 :       gel(t, j) = muluu(xa[k], xa[k+1]);
    1524     1146251 :     if (k==n) gel(t, j) = utoi(xa[k]);
    1525             :   } else {
    1526     2655025 :     for (j=1, k=1; k<n; j++, k+=2)
    1527     1372995 :       gel(t, j) = mulii(gel(xa,k), gel(xa,k+1));
    1528     1282030 :     if (k==n) gel(t, j) = icopy(gel(xa,k));
    1529             :   }
    1530     2428281 :   gel(T,1) = t;
    1531     3850996 :   for (i=2; i<=m; i++)
    1532             :   {
    1533     1422717 :     GEN u = gel(T, i-1);
    1534     1422717 :     long n = lg(u)-1;
    1535     1422717 :     t = cgetg(((n+1)>>1)+1, t_VEC);
    1536     3193338 :     for (j=1, k=1; k<n; j++, k+=2)
    1537     1770623 :       gel(t, j) = mulii(gel(u, k), gel(u, k+1));
    1538     1422715 :     if (k==n) gel(t, j) = gel(u, k);
    1539     1422715 :     gel(T, i) = t;
    1540             :   }
    1541     2428279 :   return T;
    1542             : }
    1543             : 
    1544             : /* return [A mod P[i], i=1..#P], T = ZV_producttree(P) */
    1545             : GEN
    1546    58125158 : Z_ZV_mod_tree(GEN A, GEN P, GEN T)
    1547             : {
    1548             :   long i,j,k;
    1549    58125158 :   long m = lg(T)-1, n = lg(P)-1;
    1550             :   GEN t;
    1551    58125158 :   GEN Tp = cgetg(m+1, t_VEC);
    1552    58074821 :   gel(Tp, m) = mkvec(modii(A, gmael(T,m,1)));
    1553   122782966 :   for (i=m-1; i>=1; i--)
    1554             :   {
    1555    64767933 :     GEN u = gel(T, i);
    1556    64767933 :     GEN v = gel(Tp, i+1);
    1557    64767933 :     long n = lg(u)-1;
    1558    64767933 :     t = cgetg(n+1, t_VEC);
    1559   159302574 :     for (j=1, k=1; k<n; j++, k+=2)
    1560             :     {
    1561    94618839 :       gel(t, k)   = modii(gel(v, j), gel(u, k));
    1562    94637064 :       gel(t, k+1) = modii(gel(v, j), gel(u, k+1));
    1563             :     }
    1564    64683735 :     if (k==n) gel(t, k) = gel(v, j);
    1565    64683735 :     gel(Tp, i) = t;
    1566             :   }
    1567             :   {
    1568    58015033 :     GEN u = gel(T, i+1);
    1569    58015033 :     GEN v = gel(Tp, i+1);
    1570    58015033 :     long l = lg(u)-1;
    1571    58015033 :     if (typ(P)==t_VECSMALL)
    1572             :     {
    1573    55591964 :       GEN R = cgetg(n+1, t_VECSMALL);
    1574   205574240 :       for (j=1, k=1; j<=l; j++, k+=2)
    1575             :       {
    1576   149609596 :         uel(R,k) = umodiu(gel(v, j), P[k]);
    1577   149967437 :         if (k < n)
    1578   118723305 :           uel(R,k+1) = umodiu(gel(v, j), P[k+1]);
    1579             :       }
    1580    55964644 :       return R;
    1581             :     }
    1582             :     else
    1583             :     {
    1584     2423069 :       GEN R = cgetg(n+1, t_VEC);
    1585     6623508 :       for (j=1, k=1; j<=l; j++, k+=2)
    1586             :       {
    1587     4195581 :         gel(R,k) = modii(gel(v, j), gel(P,k));
    1588     4195588 :         if (k < n)
    1589     3421913 :           gel(R,k+1) = modii(gel(v, j), gel(P,k+1));
    1590             :       }
    1591     2427927 :       return R;
    1592             :     }
    1593             :   }
    1594             : }
    1595             : 
    1596             : /* T = ZV_producttree(P), R = ZV_chinesetree(P,T) */
    1597             : GEN
    1598    39451402 : ZV_chinese_tree(GEN A, GEN P, GEN T, GEN R)
    1599             : {
    1600    39451402 :   long m = lg(T)-1, n = lg(A)-1;
    1601             :   long i,j,k;
    1602    39451402 :   GEN Tp = cgetg(m+1, t_VEC);
    1603    39443169 :   GEN M = gel(T, 1);
    1604    39443169 :   GEN t = cgetg(lg(M), t_VEC);
    1605    39395648 :   if (typ(P)==t_VECSMALL)
    1606             :   {
    1607    88464984 :     for (j=1, k=1; k<n; j++, k+=2)
    1608             :     {
    1609    64998232 :       pari_sp av = avma;
    1610    64998232 :       GEN a = mului(A[k], gel(R,k)), b = mului(A[k+1], gel(R,k+1));
    1611    64848200 :       GEN tj = modii(addii(mului(P[k],b), mului(P[k+1],a)), gel(M,j));
    1612    64956199 :       gel(t, j) = gerepileuptoint(av, tj);
    1613             :     }
    1614    23466752 :     if (k==n) gel(t, j) = modii(mului(A[k], gel(R,k)), gel(M, j));
    1615             :   } else
    1616             :   {
    1617    34015713 :     for (j=1, k=1; k<n; j++, k+=2)
    1618             :     {
    1619    18040283 :       pari_sp av = avma;
    1620    18040283 :       GEN a = mulii(gel(A,k), gel(R,k)), b = mulii(gel(A,k+1), gel(R,k+1));
    1621    18040747 :       GEN tj = modii(addii(mulii(gel(P,k),b), mulii(gel(P,k+1),a)), gel(M,j));
    1622    18073080 :       gel(t, j) = gerepileuptoint(av, tj);
    1623             :     }
    1624    15975430 :     if (k==n) gel(t, j) = modii(mulii(gel(A,k), gel(R,k)), gel(M, j));
    1625             :   }
    1626    39435502 :   gel(Tp, 1) = t;
    1627    74922360 :   for (i=2; i<=m; i++)
    1628             :   {
    1629    35453703 :     GEN u = gel(T, i-1), M = gel(T, i);
    1630    35453703 :     GEN t = cgetg(lg(M), t_VEC);
    1631    35442355 :     GEN v = gel(Tp, i-1);
    1632    35442355 :     long n = lg(v)-1;
    1633    95222162 :     for (j=1, k=1; k<n; j++, k+=2)
    1634             :     {
    1635    59735304 :       pari_sp av = avma;
    1636    59666982 :       gel(t, j) = gerepileuptoint(av, modii(addii(mulii(gel(u, k), gel(v, k+1)),
    1637    59735304 :             mulii(gel(u, k+1), gel(v, k))), gel(M, j)));
    1638             :     }
    1639    35486858 :     if (k==n) gel(t, j) = gel(v, k);
    1640    35486858 :     gel(Tp, i) = t;
    1641             :   }
    1642    39468657 :   return gmael(Tp,m,1);
    1643             : }
    1644             : 
    1645             : static GEN
    1646     1464308 : ncV_polint_center_tree(GEN vA, GEN P, GEN T, GEN R, GEN m2)
    1647             : {
    1648     1464308 :   long i, l = lg(gel(vA,1)), n = lg(P);
    1649     1464308 :   GEN mod = gmael(T, lg(T)-1, 1), V = cgetg(l, t_COL);
    1650    34972644 :   for (i=1; i < l; i++)
    1651             :   {
    1652    33508591 :     pari_sp av = avma;
    1653    33508591 :     GEN c, A = cgetg(n, typ(P));
    1654             :     long j;
    1655   199260116 :     for (j=1; j < n; j++) A[j] = mael(vA,j,i);
    1656    33466776 :     c = Fp_center(ZV_chinese_tree(A, P, T, R), mod, m2);
    1657    33492698 :     gel(V,i) = gerepileuptoint(av, c);
    1658             :   }
    1659     1464053 :   return V;
    1660             : }
    1661             : 
    1662             : static GEN
    1663      633126 : nxV_polint_center_tree(GEN vA, GEN P, GEN T, GEN R, GEN m2)
    1664             : {
    1665      633126 :   long i, j, l, n = lg(P);
    1666      633126 :   GEN mod = gmael(T, lg(T)-1, 1), V, w;
    1667      633126 :   w = cgetg(n, t_VECSMALL);
    1668     2172813 :   for(j=1; j<n; j++) w[j] = lg(gel(vA,j));
    1669      633127 :   l = vecsmall_max(w);
    1670      633121 :   V = cgetg(l, t_POL);
    1671      633043 :   V[1] = mael(vA,1,1);
    1672     3868965 :   for (i=2; i < l; i++)
    1673             :   {
    1674     3235867 :     pari_sp av = avma;
    1675     3235867 :     GEN c, A = cgetg(n, typ(P));
    1676     3235270 :     if (typ(P)==t_VECSMALL)
    1677     7769903 :       for (j=1; j < n; j++) A[j] = i < w[j] ? mael(vA,j,i): 0;
    1678             :     else
    1679     4346796 :       for (j=1; j < n; j++) gel(A,j) = i < w[j] ? gmael(vA,j,i): gen_0;
    1680     3235270 :     c = Fp_center(ZV_chinese_tree(A, P, T, R), mod, m2);
    1681     3235808 :     gel(V,i) = gerepileuptoint(av, c);
    1682             :   }
    1683      633098 :   return ZX_renormalize(V, l);
    1684             : }
    1685             : 
    1686             : static GEN
    1687        4614 : nxCV_polint_center_tree(GEN vA, GEN P, GEN T, GEN R, GEN m2)
    1688             : {
    1689        4614 :   long i, j, l = lg(gel(vA,1)), n = lg(P);
    1690        4614 :   GEN A = cgetg(n, t_VEC);
    1691        4614 :   GEN V = cgetg(l, t_COL);
    1692       90748 :   for (i=1; i < l; i++)
    1693             :   {
    1694      334652 :     for (j=1; j < n; j++) gel(A,j) = gmael(vA,j,i);
    1695       86134 :     gel(V,i) = nxV_polint_center_tree(A, P, T, R, m2);
    1696             :   }
    1697        4614 :   return V;
    1698             : }
    1699             : 
    1700             : static GEN
    1701      367742 : polint_chinese(GEN worker, GEN mA, GEN P)
    1702             : {
    1703      367742 :   long cnt, pending, n, i, j, l = lg(gel(mA,1));
    1704             :   struct pari_mt pt;
    1705             :   GEN done, va, M, A;
    1706             :   pari_timer ti;
    1707             : 
    1708      367742 :   if (l == 1) return cgetg(1, t_MAT);
    1709      338690 :   cnt = pending = 0;
    1710      338690 :   n = lg(P);
    1711      338690 :   A = cgetg(n, t_VEC);
    1712      338690 :   va = mkvec(A);
    1713      338690 :   M = cgetg(l, t_MAT);
    1714      338690 :   if (DEBUGLEVEL>4) timer_start(&ti);
    1715      338690 :   if (DEBUGLEVEL>5) err_printf("Start parallel Chinese remainder: ");
    1716      338690 :   mt_queue_start_lim(&pt, worker, l-1);
    1717     1314392 :   for (i=1; i<l || pending; i++)
    1718             :   {
    1719             :     long workid;
    1720     3692546 :     for(j=1; j < n; j++) gel(A,j) = gmael(mA,j,i);
    1721      975702 :     mt_queue_submit(&pt, i, i<l? va: NULL);
    1722      975702 :     done = mt_queue_get(&pt, &workid, &pending);
    1723      975702 :     if (done)
    1724             :     {
    1725      935750 :       gel(M,workid) = done;
    1726      935750 :       if (DEBUGLEVEL>5) err_printf("%ld%% ",(++cnt)*100/(l-1));
    1727             :     }
    1728             :   }
    1729      338690 :   if (DEBUGLEVEL>5) err_printf("\n");
    1730      338690 :   if (DEBUGLEVEL>4) timer_printf(&ti, "nmV_chinese_center");
    1731      338690 :   mt_queue_end(&pt);
    1732      338690 :   return M;
    1733             : }
    1734             : 
    1735             : GEN
    1736         840 : nxMV_polint_center_tree_worker(GEN vA, GEN T, GEN R, GEN P, GEN m2)
    1737             : {
    1738         840 :   return nxCV_polint_center_tree(vA, P, T, R, m2);
    1739             : }
    1740             : 
    1741             : static GEN
    1742         430 : nxMV_polint_center_tree_seq(GEN vA, GEN P, GEN T, GEN R, GEN m2)
    1743             : {
    1744         430 :   long i, j, l = lg(gel(vA,1)), n = lg(P);
    1745         430 :   GEN A = cgetg(n, t_VEC);
    1746         430 :   GEN V = cgetg(l, t_MAT);
    1747        4204 :   for (i=1; i < l; i++)
    1748             :   {
    1749       15299 :     for (j=1; j < n; j++) gel(A,j) = gmael(vA,j,i);
    1750        3774 :     gel(V,i) = nxCV_polint_center_tree(A, P, T, R, m2);
    1751             :   }
    1752         430 :   return V;
    1753             : }
    1754             : 
    1755             : static GEN
    1756          90 : nxMV_polint_center_tree(GEN mA, GEN P, GEN T, GEN R, GEN m2)
    1757             : {
    1758          90 :   GEN worker = snm_closure(is_entry("_nxMV_polint_worker"), mkvec4(T, R, P, m2));
    1759          90 :   return polint_chinese(worker, mA, P);
    1760             : }
    1761             : 
    1762             : static GEN
    1763      111805 : nmV_polint_center_tree_seq(GEN vA, GEN P, GEN T, GEN R, GEN m2)
    1764             : {
    1765      111805 :   long i, j, l = lg(gel(vA,1)), n = lg(P);
    1766      111805 :   GEN A = cgetg(n, t_VEC);
    1767      111806 :   GEN V = cgetg(l, t_MAT);
    1768      625729 :   for (i=1; i < l; i++)
    1769             :   {
    1770     3001010 :     for (j=1; j < n; j++) gel(A,j) = gmael(vA,j,i);
    1771      513926 :     gel(V,i) = ncV_polint_center_tree(A, P, T, R, m2);
    1772             :   }
    1773      111803 :   return V;
    1774             : }
    1775             : 
    1776             : GEN
    1777      934764 : nmV_polint_center_tree_worker(GEN vA, GEN T, GEN R, GEN P, GEN m2)
    1778             : {
    1779      934764 :   return ncV_polint_center_tree(vA, P, T, R, m2);
    1780             : }
    1781             : 
    1782             : static GEN
    1783      367651 : nmV_polint_center_tree(GEN mA, GEN P, GEN T, GEN R, GEN m2)
    1784             : {
    1785      367651 :   GEN worker = snm_closure(is_entry("_polint_worker"), mkvec4(T, R, P, m2));
    1786      367652 :   return polint_chinese(worker, mA, P);
    1787             : }
    1788             : 
    1789             : /* return [A mod P[i], i=1..#P] */
    1790             : GEN
    1791           0 : Z_ZV_mod(GEN A, GEN P)
    1792             : {
    1793           0 :   pari_sp av = avma;
    1794           0 :   return gerepilecopy(av, Z_ZV_mod_tree(A, P, ZV_producttree(P)));
    1795             : }
    1796             : /* P a t_VECSMALL */
    1797             : GEN
    1798           0 : Z_nv_mod(GEN A, GEN P)
    1799             : {
    1800           0 :   pari_sp av = avma;
    1801           0 :   return gerepileuptoleaf(av, Z_ZV_mod_tree(A, P, ZV_producttree(P)));
    1802             : }
    1803             : /* B a ZX, T = ZV_producttree(P) */
    1804             : GEN
    1805     1462854 : ZX_nv_mod_tree(GEN B, GEN A, GEN T)
    1806             : {
    1807             :   pari_sp av;
    1808     1462854 :   long i, j, l = lg(B), n = lg(A)-1;
    1809     1462854 :   GEN V = cgetg(n+1, t_VEC);
    1810     7168747 :   for (j=1; j <= n; j++)
    1811             :   {
    1812     5706035 :     gel(V, j) = cgetg(l, t_VECSMALL);
    1813     5705931 :     mael(V, j, 1) = B[1]&VARNBITS;
    1814             :   }
    1815     1462712 :   av = avma;
    1816    12775237 :   for (i=2; i < l; i++)
    1817             :   {
    1818    11313111 :     GEN v = Z_ZV_mod_tree(gel(B, i), A, T);
    1819    78684568 :     for (j=1; j <= n; j++)
    1820    67377326 :       mael(V, j, i) = v[j];
    1821    11307242 :     set_avma(av);
    1822             :   }
    1823     7168495 :   for (j=1; j <= n; j++)
    1824     5706334 :     (void) Flx_renormalize(gel(V, j), l);
    1825     1462161 :   return V;
    1826             : }
    1827             : 
    1828             : static GEN
    1829      299615 : to_ZX(GEN a, long v) { return typ(a)==t_INT? scalarpol(a,v): a; }
    1830             : 
    1831             : GEN
    1832       31091 : ZXX_nv_mod_tree(GEN P, GEN xa, GEN T, long w)
    1833             : {
    1834       31091 :   pari_sp av = avma;
    1835       31091 :   long i, j, l = lg(P), n = lg(xa)-1;
    1836       31091 :   GEN V = cgetg(n+1, t_VEC);
    1837      126473 :   for (j=1; j <= n; j++)
    1838             :   {
    1839       95382 :     gel(V, j) = cgetg(l, t_POL);
    1840       95382 :     mael(V, j, 1) = P[1]&VARNBITS;
    1841             :   }
    1842      249467 :   for (i=2; i < l; i++)
    1843             :   {
    1844      218376 :     GEN v = ZX_nv_mod_tree(to_ZX(gel(P, i), w), xa, T);
    1845      868278 :     for (j=1; j <= n; j++)
    1846      649902 :       gmael(V, j, i) = gel(v,j);
    1847             :   }
    1848      126473 :   for (j=1; j <= n; j++)
    1849       95382 :     (void) FlxX_renormalize(gel(V, j), l);
    1850       31091 :   return gerepilecopy(av, V);
    1851             : }
    1852             : 
    1853             : GEN
    1854        4048 : ZXC_nv_mod_tree(GEN C, GEN xa, GEN T, long w)
    1855             : {
    1856        4048 :   pari_sp av = avma;
    1857        4048 :   long i, j, l = lg(C), n = lg(xa)-1;
    1858        4048 :   GEN V = cgetg(n+1, t_VEC);
    1859       16952 :   for (j = 1; j <= n; j++)
    1860       12904 :     gel(V, j) = cgetg(l, t_COL);
    1861       85270 :   for (i = 1; i < l; i++)
    1862             :   {
    1863       81229 :     GEN v = ZX_nv_mod_tree(to_ZX(gel(C, i), w), xa, T);
    1864      359594 :     for (j = 1; j <= n; j++)
    1865      278372 :       gmael(V, j, i) = gel(v,j);
    1866             :   }
    1867        4041 :   return gerepilecopy(av, V);
    1868             : }
    1869             : 
    1870             : GEN
    1871         430 : ZXM_nv_mod_tree(GEN M, GEN xa, GEN T, long w)
    1872             : {
    1873         430 :   pari_sp av = avma;
    1874         430 :   long i, j, l = lg(M), n = lg(xa)-1;
    1875         430 :   GEN V = cgetg(n+1, t_VEC);
    1876        2083 :   for (j=1; j <= n; j++)
    1877        1653 :     gel(V, j) = cgetg(l, t_MAT);
    1878        4204 :   for (i=1; i < l; i++)
    1879             :   {
    1880        3774 :     GEN v = ZXC_nv_mod_tree(gel(M, i), xa, T, w);
    1881       15299 :     for (j=1; j <= n; j++)
    1882       11525 :       gmael(V, j, i) = gel(v,j);
    1883             :   }
    1884         430 :   return gerepilecopy(av, V);
    1885             : }
    1886             : 
    1887             : GEN
    1888     1208396 : ZV_nv_mod_tree(GEN B, GEN A, GEN T)
    1889             : {
    1890             :   pari_sp av;
    1891     1208396 :   long i, j, l = lg(B), n = lg(A)-1;
    1892     1208396 :   GEN V = cgetg(n+1, t_VEC);
    1893     6207533 :   for (j=1; j <= n; j++) gel(V, j) = cgetg(l, t_VECSMALL);
    1894     1208299 :   av = avma;
    1895    45584389 :   for (i=1; i < l; i++)
    1896             :   {
    1897    44379500 :     GEN v = Z_ZV_mod_tree(gel(B, i), A, T);
    1898   246877933 :     for (j=1; j <= n; j++) mael(V, j, i) = v[j];
    1899    44352435 :     set_avma(av);
    1900             :   }
    1901     1204889 :   return V;
    1902             : }
    1903             : 
    1904             : GEN
    1905      178783 : ZM_nv_mod_tree(GEN M, GEN xa, GEN T)
    1906             : {
    1907      178783 :   pari_sp av = avma;
    1908      178783 :   long i, j, l = lg(M), n = lg(xa)-1;
    1909      178783 :   GEN V = cgetg(n+1, t_VEC);
    1910      894577 :   for (j=1; j <= n; j++) gel(V, j) = cgetg(l, t_MAT);
    1911     1386949 :   for (i=1; i < l; i++)
    1912             :   {
    1913     1208165 :     GEN v = ZV_nv_mod_tree(gel(M, i), xa, T);
    1914     6207087 :     for (j=1; j <= n; j++) gmael(V, j, i) = gel(v,j);
    1915             :   }
    1916      178784 :   return gerepilecopy(av, V);
    1917             : }
    1918             : 
    1919             : static GEN
    1920     2424072 : ZV_sqr(GEN z)
    1921             : {
    1922     2424072 :   long i,l = lg(z);
    1923     2424072 :   GEN x = cgetg(l, t_VEC);
    1924     2424071 :   if (typ(z)==t_VECSMALL)
    1925     5659296 :     for (i=1; i<l; i++) gel(x,i) = sqru(z[i]);
    1926             :   else
    1927     4356610 :     for (i=1; i<l; i++) gel(x,i) = sqri(gel(z,i));
    1928     2424074 :   return x;
    1929             : }
    1930             : 
    1931             : static GEN
    1932    12440874 : ZT_sqr(GEN x)
    1933             : {
    1934    12440874 :   if (typ(x) == t_INT) return sqri(x);
    1935    16281234 :   pari_APPLY_type(t_VEC, ZT_sqr(gel(x,i)))
    1936             : }
    1937             : 
    1938             : static GEN
    1939     2424073 : ZV_invdivexact(GEN y, GEN x)
    1940             : {
    1941     2424073 :   long i, l = lg(y);
    1942     2424073 :   GEN z = cgetg(l,t_VEC);
    1943     2424071 :   if (typ(x)==t_VECSMALL)
    1944     5659144 :     for (i=1; i<l; i++)
    1945             :     {
    1946     4513247 :       pari_sp av = avma;
    1947     4513247 :       ulong a = Fl_inv(umodiu(diviuexact(gel(y,i),x[i]), x[i]), x[i]);
    1948     4513419 :       set_avma(av); gel(z,i) = utoi(a);
    1949             :     }
    1950             :   else
    1951     4356609 :     for (i=1; i<l; i++)
    1952     3078440 :       gel(z,i) = Fp_inv(diviiexact(gel(y,i), gel(x,i)), gel(x,i));
    1953     2424066 :   return z;
    1954             : }
    1955             : 
    1956             : /* P t_VECSMALL or t_VEC of t_INT  */
    1957             : GEN
    1958     2424067 : ZV_chinesetree(GEN P, GEN T)
    1959             : {
    1960     2424067 :   GEN T2 = ZT_sqr(T), P2 = ZV_sqr(P);
    1961     2424068 :   GEN mod = gmael(T,lg(T)-1,1);
    1962     2424068 :   return ZV_invdivexact(Z_ZV_mod_tree(mod, P2, T2), P);
    1963             : }
    1964             : 
    1965             : static GEN
    1966      928730 : gc_chinese(pari_sp av, GEN T, GEN a, GEN *pt_mod)
    1967             : {
    1968      928730 :   if (!pt_mod)
    1969       12237 :     return gerepileupto(av, a);
    1970             :   else
    1971             :   {
    1972      916493 :     GEN mod = gmael(T, lg(T)-1, 1);
    1973      916493 :     gerepileall(av, 2, &a, &mod);
    1974      916494 :     *pt_mod = mod;
    1975      916494 :     return a;
    1976             :   }
    1977             : }
    1978             : 
    1979             : GEN
    1980      157964 : ZV_chinese_center(GEN A, GEN P, GEN *pt_mod)
    1981             : {
    1982      157964 :   pari_sp av = avma;
    1983      157964 :   GEN T = ZV_producttree(P);
    1984      157964 :   GEN R = ZV_chinesetree(P, T);
    1985      157964 :   GEN a = ZV_chinese_tree(A, P, T, R);
    1986      157964 :   GEN mod = gmael(T, lg(T)-1, 1);
    1987      157964 :   GEN ca = Fp_center(a, mod, shifti(mod,-1));
    1988      157964 :   return gc_chinese(av, T, ca, pt_mod);
    1989             : }
    1990             : 
    1991             : GEN
    1992        5110 : ZV_chinese(GEN A, GEN P, GEN *pt_mod)
    1993             : {
    1994        5110 :   pari_sp av = avma;
    1995        5110 :   GEN T = ZV_producttree(P);
    1996        5110 :   GEN R = ZV_chinesetree(P, T);
    1997        5110 :   GEN a = ZV_chinese_tree(A, P, T, R);
    1998        5110 :   return gc_chinese(av, T, a, pt_mod);
    1999             : }
    2000             : 
    2001             : GEN
    2002      159226 : nxV_chinese_center_tree(GEN A, GEN P, GEN T, GEN R)
    2003             : {
    2004      159226 :   pari_sp av = avma;
    2005      159226 :   GEN m2 = shifti(gmael(T, lg(T)-1, 1), -1);
    2006      159226 :   GEN a = nxV_polint_center_tree(A, P, T, R, m2);
    2007      159227 :   return gerepileupto(av, a);
    2008             : }
    2009             : 
    2010             : GEN
    2011      387722 : nxV_chinese_center(GEN A, GEN P, GEN *pt_mod)
    2012             : {
    2013      387722 :   pari_sp av = avma;
    2014      387722 :   GEN T = ZV_producttree(P);
    2015      387722 :   GEN R = ZV_chinesetree(P, T);
    2016      387720 :   GEN m2 = shifti(gmael(T, lg(T)-1, 1), -1);
    2017      387720 :   GEN a = nxV_polint_center_tree(A, P, T, R, m2);
    2018      387721 :   return gc_chinese(av, T, a, pt_mod);
    2019             : }
    2020             : 
    2021             : GEN
    2022       10193 : ncV_chinese_center(GEN A, GEN P, GEN *pt_mod)
    2023             : {
    2024       10193 :   pari_sp av = avma;
    2025       10193 :   GEN T = ZV_producttree(P);
    2026       10193 :   GEN R = ZV_chinesetree(P, T);
    2027       10193 :   GEN m2 = shifti(gmael(T, lg(T)-1, 1), -1);
    2028       10193 :   GEN a = ncV_polint_center_tree(A, P, T, R, m2);
    2029       10193 :   return gc_chinese(av, T, a, pt_mod);
    2030             : }
    2031             : 
    2032             : GEN
    2033        5457 : ncV_chinese_center_tree(GEN A, GEN P, GEN T, GEN R)
    2034             : {
    2035        5457 :   pari_sp av = avma;
    2036        5457 :   GEN m2 = shifti(gmael(T, lg(T)-1, 1), -1);
    2037        5457 :   GEN a = ncV_polint_center_tree(A, P, T, R, m2);
    2038        5457 :   return gerepileupto(av, a);
    2039             : }
    2040             : 
    2041             : GEN
    2042           0 : nmV_chinese_center_tree(GEN A, GEN P, GEN T, GEN R)
    2043             : {
    2044           0 :   pari_sp av = avma;
    2045           0 :   GEN m2 = shifti(gmael(T, lg(T)-1, 1), -1);
    2046           0 :   GEN a = nmV_polint_center_tree(A, P, T, R, m2);
    2047           0 :   return gerepileupto(av, a);
    2048             : }
    2049             : 
    2050             : GEN
    2051      111806 : nmV_chinese_center_tree_seq(GEN A, GEN P, GEN T, GEN R)
    2052             : {
    2053      111806 :   pari_sp av = avma;
    2054      111806 :   GEN m2 = shifti(gmael(T, lg(T)-1, 1), -1);
    2055      111805 :   GEN a = nmV_polint_center_tree_seq(A, P, T, R, m2);
    2056      111803 :   return gerepileupto(av, a);
    2057             : }
    2058             : 
    2059             : GEN
    2060      367652 : nmV_chinese_center(GEN A, GEN P, GEN *pt_mod)
    2061             : {
    2062      367652 :   pari_sp av = avma;
    2063      367652 :   GEN T = ZV_producttree(P);
    2064      367652 :   GEN R = ZV_chinesetree(P, T);
    2065      367651 :   GEN m2 = shifti(gmael(T, lg(T)-1, 1), -1);
    2066      367651 :   GEN a = nmV_polint_center_tree(A, P, T, R, m2);
    2067      367652 :   return gc_chinese(av, T, a, pt_mod);
    2068             : }
    2069             : 
    2070             : GEN
    2071           0 : nxCV_chinese_center_tree(GEN A, GEN P, GEN T, GEN R)
    2072             : {
    2073           0 :   pari_sp av = avma;
    2074           0 :   GEN m2 = shifti(gmael(T, lg(T)-1, 1), -1);
    2075           0 :   GEN a = nxCV_polint_center_tree(A, P, T, R, m2);
    2076           0 :   return gerepileupto(av, a);
    2077             : }
    2078             : 
    2079             : GEN
    2080           0 : nxCV_chinese_center(GEN A, GEN P, GEN *pt_mod)
    2081             : {
    2082           0 :   pari_sp av = avma;
    2083           0 :   GEN T = ZV_producttree(P);
    2084           0 :   GEN R = ZV_chinesetree(P, T);
    2085           0 :   GEN m2 = shifti(gmael(T, lg(T)-1, 1), -1);
    2086           0 :   GEN a = nxCV_polint_center_tree(A, P, T, R, m2);
    2087           0 :   return gc_chinese(av, T, a, pt_mod);
    2088             : }
    2089             : 
    2090             : GEN
    2091         430 : nxMV_chinese_center_tree_seq(GEN A, GEN P, GEN T, GEN R)
    2092             : {
    2093         430 :   pari_sp av = avma;
    2094         430 :   GEN m2 = shifti(gmael(T, lg(T)-1, 1), -1);
    2095         430 :   GEN a = nxMV_polint_center_tree_seq(A, P, T, R, m2);
    2096         430 :   return gerepileupto(av, a);
    2097             : }
    2098             : 
    2099             : GEN
    2100          90 : nxMV_chinese_center(GEN A, GEN P, GEN *pt_mod)
    2101             : {
    2102          90 :   pari_sp av = avma;
    2103          90 :   GEN T = ZV_producttree(P);
    2104          90 :   GEN R = ZV_chinesetree(P, T);
    2105          90 :   GEN m2 = shifti(gmael(T, lg(T)-1, 1), -1);
    2106          90 :   GEN a = nxMV_polint_center_tree(A, P, T, R, m2);
    2107          90 :   return gc_chinese(av, T, a, pt_mod);
    2108             : }
    2109             : 
    2110             : /**********************************************************************
    2111             :  **                    Powering  over (Z/NZ)^*, small N              **
    2112             :  **********************************************************************/
    2113             : 
    2114             : /* 2^n mod p; assume n > 1 */
    2115             : static ulong
    2116    11955777 : Fl_2powu_pre(ulong n, ulong p, ulong pi)
    2117             : {
    2118    11955777 :   ulong y = 2;
    2119    11955777 :   int j = 1+bfffo(n);
    2120             :   /* normalize, i.e set highest bit to 1 (we know n != 0) */
    2121    11955777 :   n<<=j; j = BITS_IN_LONG-j; /* first bit is now implicit */
    2122   522082083 :   for (; j; n<<=1,j--)
    2123             :   {
    2124   510115228 :     y = Fl_sqr_pre(y,p,pi);
    2125   510110065 :     if (n & HIGHBIT) y = Fl_double(y, p);
    2126             :   }
    2127    11966855 :   return y;
    2128             : }
    2129             : 
    2130             : /* 2^n mod p; assume n > 1 and !(p & HIGHMASK) */
    2131             : static ulong
    2132     4160241 : Fl_2powu(ulong n, ulong p)
    2133             : {
    2134     4160241 :   ulong y = 2;
    2135     4160241 :   int j = 1+bfffo(n);
    2136             :   /* normalize, i.e set highest bit to 1 (we know n != 0) */
    2137     4160241 :   n<<=j; j = BITS_IN_LONG-j; /* first bit is now implicit */
    2138    23247680 :   for (; j; n<<=1,j--)
    2139             :   {
    2140    19087449 :     y = (y*y) % p;
    2141    19087449 :     if (n & HIGHBIT) y = Fl_double(y, p);
    2142             :   }
    2143     4160231 :   return y;
    2144             : }
    2145             : 
    2146             : /* allow pi = 0 */
    2147             : ulong
    2148    98102438 : Fl_powu_pre(ulong x, ulong n0, ulong p, ulong pi)
    2149             : {
    2150             :   ulong y, z, n;
    2151    98102438 :   if (!pi) return Fl_powu(x, n0, p);
    2152    95663569 :   if (n0 <= 1)
    2153             :   { /* frequent special cases */
    2154     6999230 :     if (n0 == 1) return x;
    2155      103726 :     if (n0 == 0) return 1;
    2156             :   }
    2157    88664331 :   if (x <= 2)
    2158             :   {
    2159    12213950 :     if (x == 2) return Fl_2powu_pre(n0, p, pi);
    2160      256970 :     return x; /* 0 or 1 */
    2161             :   }
    2162    76450381 :   y = 1; z = x; n = n0;
    2163             :   for(;;)
    2164             :   {
    2165   540423882 :     if (n&1) y = Fl_mul_pre(y,z,p,pi);
    2166   540790848 :     n>>=1; if (!n) return y;
    2167   464155502 :     z = Fl_sqr_pre(z,p,pi);
    2168             :   }
    2169             : }
    2170             : 
    2171             : ulong
    2172   146289334 : Fl_powu(ulong x, ulong n0, ulong p)
    2173             : {
    2174             :   ulong y, z, n;
    2175   146289334 :   if (n0 <= 2)
    2176             :   { /* frequent special cases */
    2177    65647154 :     if (n0 == 2) return Fl_sqr(x,p);
    2178    32108389 :     if (n0 == 1) return x;
    2179     1911345 :     if (n0 == 0) return 1;
    2180             :   }
    2181    80614087 :   if (x <= 1) return x; /* 0 or 1 */
    2182    80182470 :   if (p & HIGHMASK)
    2183     7826091 :     return Fl_powu_pre(x, n0, p, get_Fl_red(p));
    2184    72356379 :   if (x == 2) return Fl_2powu(n0, p);
    2185    68196140 :   y = 1; z = x; n = n0;
    2186             :   for(;;)
    2187             :   {
    2188   370067861 :     if (n&1) y = (y*z) % p;
    2189   370067861 :     n>>=1; if (!n) return y;
    2190   301871721 :     z = (z*z) % p;
    2191             :   }
    2192             : }
    2193             : 
    2194             : /* Reduce data dependency to maximize internal parallelism; allow pi = 0 */
    2195             : GEN
    2196    12769114 : Fl_powers_pre(ulong x, long n, ulong p, ulong pi)
    2197             : {
    2198             :   long i, k;
    2199    12769114 :   GEN z = cgetg(n + 2, t_VECSMALL);
    2200    12761483 :   z[1] = 1; if (n == 0) return z;
    2201    12761483 :   z[2] = x;
    2202    12761483 :   if (pi)
    2203             :   {
    2204    90032826 :     for (i = 3, k=2; i <= n; i+=2, k++)
    2205             :     {
    2206    77480559 :       z[i] = Fl_sqr_pre(z[k], p, pi);
    2207    77487942 :       z[i+1] = Fl_mul_pre(z[k], z[k+1], p, pi);
    2208             :     }
    2209    12552267 :     if (i==n+1) z[i] = Fl_sqr_pre(z[k], p, pi);
    2210             :   }
    2211      212770 :   else if (p & HIGHMASK)
    2212             :   {
    2213           0 :     for (i = 3, k=2; i <= n; i+=2, k++)
    2214             :     {
    2215           0 :       z[i] = Fl_sqr(z[k], p);
    2216           0 :       z[i+1] = Fl_mul(z[k], z[k+1], p);
    2217             :     }
    2218           0 :     if (i==n+1) z[i] = Fl_sqr(z[k], p);
    2219             :   }
    2220             :   else
    2221   400516205 :     for (i = 2; i <= n; i++) z[i+1] = (z[i] * x) % p;
    2222    12767505 :   return z;
    2223             : }
    2224             : 
    2225             : GEN
    2226      295917 : Fl_powers(ulong x, long n, ulong p)
    2227             : {
    2228      295917 :   return Fl_powers_pre(x, n, p, (p & HIGHMASK)? get_Fl_red(p): 0);
    2229             : }
    2230             : 
    2231             : /**********************************************************************
    2232             :  **                    Powering  over (Z/NZ)^*, large N              **
    2233             :  **********************************************************************/
    2234             : typedef struct muldata {
    2235             :   GEN (*sqr)(void * E, GEN x);
    2236             :   GEN (*mul)(void * E, GEN x, GEN y);
    2237             :   GEN (*mul2)(void * E, GEN x);
    2238             : } muldata;
    2239             : 
    2240             : /* modified Barrett reduction with one fold */
    2241             : /* See Fast Modular Reduction, W. Hasenplaugh, G. Gaubatz, V. Gopal, ARITH 18 */
    2242             : 
    2243             : static GEN
    2244       15011 : Fp_invmBarrett(GEN p, long s)
    2245             : {
    2246       15011 :   GEN R, Q = dvmdii(int2n(3*s),p,&R);
    2247       15011 :   return mkvec2(Q,R);
    2248             : }
    2249             : 
    2250             : /* a <= (N-1)^2, 2^(2s-2) <= N < 2^(2s). Return 0 <= r < N such that
    2251             :  * a = r (mod N) */
    2252             : static GEN
    2253     9222784 : Fp_rem_mBarrett(GEN a, GEN B, long s, GEN N)
    2254             : {
    2255     9222784 :   pari_sp av = avma;
    2256     9222784 :   GEN P = gel(B, 1), Q = gel(B, 2); /* 2^(3s) = P N + Q, 0 <= Q < N */
    2257     9222784 :   long t = expi(P)+1; /* 2^(t-1) <= P < 2^t */
    2258     9222784 :   GEN u = shifti(a, -3*s), v = remi2n(a, 3*s); /* a = 2^(3s)u + v */
    2259     9222784 :   GEN A = addii(v, mulii(Q,u)); /* 0 <= A < 2^(3s+1) */
    2260     9222784 :   GEN q = shifti(mulii(shifti(A, t-3*s), P), -t); /* A/N - 4 < q <= A/N */
    2261     9222784 :   GEN r = subii(A, mulii(q, N));
    2262     9222784 :   GEN sr= subii(r,N);     /* 0 <= r < 4*N */
    2263     9222784 :   if (signe(sr)<0) return gerepileuptoint(av, r);
    2264     5036105 :   r=sr; sr = subii(r,N);  /* 0 <= r < 3*N */
    2265     5036105 :   if (signe(sr)<0) return gerepileuptoint(av, r);
    2266      187201 :   r=sr; sr = subii(r,N);  /* 0 <= r < 2*N */
    2267      187201 :   return gerepileuptoint(av, signe(sr)>=0 ? sr:r);
    2268             : }
    2269             : 
    2270             : /* Montgomery reduction */
    2271             : 
    2272             : INLINE ulong
    2273      694374 : init_montdata(GEN N) { return (ulong) -invmod2BIL(mod2BIL(N)); }
    2274             : 
    2275             : struct montred
    2276             : {
    2277             :   GEN N;
    2278             :   ulong inv;
    2279             : };
    2280             : 
    2281             : /* Montgomery reduction */
    2282             : static GEN
    2283    69080027 : _sqr_montred(void * E, GEN x)
    2284             : {
    2285    69080027 :   struct montred * D = (struct montred *) E;
    2286    69080027 :   return red_montgomery(sqri(x), D->N, D->inv);
    2287             : }
    2288             : 
    2289             : /* Montgomery reduction */
    2290             : static GEN
    2291     7108970 : _mul_montred(void * E, GEN x, GEN y)
    2292             : {
    2293     7108970 :   struct montred * D = (struct montred *) E;
    2294     7108970 :   return red_montgomery(mulii(x, y), D->N, D->inv);
    2295             : }
    2296             : 
    2297             : static GEN
    2298    11250690 : _mul2_montred(void * E, GEN x)
    2299             : {
    2300    11250690 :   struct montred * D = (struct montred *) E;
    2301    11250690 :   GEN z = shifti(_sqr_montred(E, x), 1);
    2302    11248550 :   long l = lgefint(D->N);
    2303    11886063 :   while (lgefint(z) > l) z = subii(z, D->N);
    2304    11248941 :   return z;
    2305             : }
    2306             : 
    2307             : static GEN
    2308    23162046 : _sqr_remii(void* N, GEN x)
    2309    23162046 : { return remii(sqri(x), (GEN) N); }
    2310             : 
    2311             : static GEN
    2312     1523445 : _mul_remii(void* N, GEN x, GEN y)
    2313     1523445 : { return remii(mulii(x, y), (GEN) N); }
    2314             : 
    2315             : static GEN
    2316     3210061 : _mul2_remii(void* N, GEN x)
    2317     3210061 : { return Fp_double(_sqr_remii(N, x), (GEN)N); }
    2318             : 
    2319             : struct redbarrett
    2320             : {
    2321             :   GEN iM, N;
    2322             :   long s;
    2323             : };
    2324             : 
    2325             : static GEN
    2326     8446594 : _sqr_remiibar(void *E, GEN x)
    2327             : {
    2328     8446594 :   struct redbarrett * D = (struct redbarrett *) E;
    2329     8446594 :   return Fp_rem_mBarrett(sqri(x), D->iM, D->s, D->N);
    2330             : }
    2331             : 
    2332             : static GEN
    2333      776190 : _mul_remiibar(void *E, GEN x, GEN y)
    2334             : {
    2335      776190 :   struct redbarrett * D = (struct redbarrett *) E;
    2336      776190 :   return Fp_rem_mBarrett(mulii(x, y), D->iM, D->s, D->N);
    2337             : }
    2338             : 
    2339             : static GEN
    2340     2078975 : _mul2_remiibar(void *E, GEN x)
    2341             : {
    2342     2078975 :   struct redbarrett * D = (struct redbarrett *) E;
    2343     2078975 :   return Fp_double(_sqr_remiibar(E, x), D->N);
    2344             : }
    2345             : 
    2346             : static long
    2347      886229 : Fp_select_red(GEN *y, ulong k, GEN N, long lN, muldata *D, void **pt_E)
    2348             : {
    2349      886229 :   if (lN >= Fp_POW_BARRETT_LIMIT && (k==0 || ((double)k)*expi(*y) > 2 + expi(N)))
    2350             :   {
    2351       15011 :     struct redbarrett * E = (struct redbarrett *) stack_malloc(sizeof(struct redbarrett));
    2352       15011 :     D->sqr = &_sqr_remiibar;
    2353       15011 :     D->mul = &_mul_remiibar;
    2354       15011 :     D->mul2 = &_mul2_remiibar;
    2355       15011 :     E->N = N;
    2356       15011 :     E->s = 1+(expi(N)>>1);
    2357       15011 :     E->iM = Fp_invmBarrett(N, E->s);
    2358       15011 :     *pt_E = (void*) E;
    2359       15011 :     return 0;
    2360             :   }
    2361      871218 :   else if (mod2(N) && lN < Fp_POW_REDC_LIMIT)
    2362             :   {
    2363      694367 :     struct montred * E = (struct montred *) stack_malloc(sizeof(struct montred));
    2364      694368 :     *y = remii(shifti(*y, bit_accuracy(lN)), N);
    2365      694374 :     D->sqr = &_sqr_montred;
    2366      694374 :     D->mul = &_mul_montred;
    2367      694374 :     D->mul2 = &_mul2_montred;
    2368      694374 :     E->N = N;
    2369      694374 :     E->inv = init_montdata(N);
    2370      694371 :     *pt_E = (void*) E;
    2371      694371 :     return 1;
    2372             :   }
    2373             :   else
    2374             :   {
    2375      176851 :     D->sqr = &_sqr_remii;
    2376      176851 :     D->mul = &_mul_remii;
    2377      176851 :     D->mul2 = &_mul2_remii;
    2378      176851 :     *pt_E = (void*) N;
    2379      176851 :     return 0;
    2380             :   }
    2381             : }
    2382             : 
    2383             : GEN
    2384     2868174 : Fp_powu(GEN A, ulong k, GEN N)
    2385             : {
    2386     2868174 :   long lN = lgefint(N);
    2387             :   int base_is_2, use_montgomery;
    2388             :   muldata D;
    2389             :   void *E;
    2390             :   pari_sp av;
    2391             : 
    2392     2868174 :   if (lN == 3) {
    2393     1438634 :     ulong n = uel(N,2);
    2394     1438634 :     return utoi( Fl_powu(umodiu(A, n), k, n) );
    2395             :   }
    2396     1429540 :   if (k <= 2)
    2397             :   { /* frequent special cases */
    2398      814168 :     if (k == 2) return Fp_sqr(A,N);
    2399      268328 :     if (k == 1) return A;
    2400           0 :     if (k == 0) return gen_1;
    2401             :   }
    2402      615372 :   av = avma; A = modii(A,N);
    2403      615372 :   base_is_2 = 0;
    2404      615372 :   if (lgefint(A) == 3) switch(A[2])
    2405             :   {
    2406         860 :     case 1: set_avma(av); return gen_1;
    2407       42678 :     case 2:  base_is_2 = 1; break;
    2408             :   }
    2409             : 
    2410             :   /* TODO: Move this out of here and use for general modular computations */
    2411      614512 :   use_montgomery = Fp_select_red(&A, k, N, lN, &D, &E);
    2412      614512 :   if (base_is_2)
    2413       42678 :     A = gen_powu_fold_i(A, k, E, D.sqr, D.mul2);
    2414             :   else
    2415      571834 :     A = gen_powu_i(A, k, E, D.sqr, D.mul);
    2416      614512 :   if (use_montgomery)
    2417             :   {
    2418      523639 :     A = red_montgomery(A, N, ((struct montred *) E)->inv);
    2419      523639 :     if (cmpii(A, N) >= 0) A = subii(A,N);
    2420             :   }
    2421      614512 :   return gerepileuptoint(av, A);
    2422             : }
    2423             : 
    2424             : GEN
    2425       29327 : Fp_pows(GEN A, long k, GEN N)
    2426             : {
    2427       29327 :   if (lgefint(N) == 3) {
    2428       14869 :     ulong n = N[2];
    2429       14869 :     ulong a = umodiu(A, n);
    2430       14869 :     if (k < 0) {
    2431         133 :       a = Fl_inv(a, n);
    2432         133 :       k = -k;
    2433             :     }
    2434       14869 :     return utoi( Fl_powu(a, (ulong)k, n) );
    2435             :   }
    2436       14458 :   if (k < 0) { A = Fp_inv(A, N); k = -k; };
    2437       14458 :   return Fp_powu(A, (ulong)k, N);
    2438             : }
    2439             : 
    2440             : /* A^K mod N */
    2441             : GEN
    2442    35823517 : Fp_pow(GEN A, GEN K, GEN N)
    2443             : {
    2444             :   pari_sp av;
    2445    35823517 :   long s, lN = lgefint(N), sA, sy;
    2446             :   int base_is_2, use_montgomery;
    2447             :   GEN y;
    2448             :   muldata D;
    2449             :   void *E;
    2450             : 
    2451    35823517 :   s = signe(K);
    2452    35823517 :   if (!s) return dvdii(A,N)? gen_0: gen_1;
    2453    34793458 :   if (lN == 3 && lgefint(K) == 3)
    2454             :   {
    2455    34032077 :     ulong n = N[2], a = umodiu(A, n);
    2456    34032408 :     if (s < 0) a = Fl_inv(a, n);
    2457    34032495 :     if (a <= 1) return utoi(a); /* 0 or 1 */
    2458    30518263 :     return utoi(Fl_powu(a, uel(K,2), n));
    2459             :   }
    2460             : 
    2461      761381 :   av = avma;
    2462      761381 :   if (s < 0) y = Fp_inv(A,N);
    2463             :   else
    2464             :   {
    2465      759444 :     y = modii(A,N);
    2466      759663 :     if (!signe(y)) { set_avma(av); return gen_0; }
    2467             :   }
    2468      761600 :   if (lgefint(K) == 3) return gerepileuptoint(av, Fp_powu(y, K[2], N));
    2469             : 
    2470      271919 :   base_is_2 = 0;
    2471      271919 :   sy = abscmpii(y, shifti(N,-1)) > 0;
    2472      271913 :   if (sy) y = subii(N,y);
    2473      271918 :   sA = sy && mod2(K);
    2474      271917 :   if (lgefint(y) == 3) switch(y[2])
    2475             :   {
    2476         207 :     case 1:  set_avma(av); return sA ? subis(N,1): gen_1;
    2477      147140 :     case 2:  base_is_2 = 1; break;
    2478             :   }
    2479             : 
    2480             :   /* TODO: Move this out of here and use for general modular computations */
    2481      271710 :   use_montgomery = Fp_select_red(&y, 0UL, N, lN, &D, &E);
    2482      271722 :   if (base_is_2)
    2483      147153 :     y = gen_pow_fold_i(y, K, E, D.sqr, D.mul2);
    2484             :   else
    2485      124569 :     y = gen_pow_i(y, K, E, D.sqr, D.mul);
    2486      271736 :   if (use_montgomery)
    2487             :   {
    2488      170738 :     y = red_montgomery(y, N, ((struct montred *) E)->inv);
    2489      170740 :     if (cmpii(y,N) >= 0) y = subii(y,N);
    2490             :   }
    2491      271738 :   if (sA) y = subii(N, y);
    2492      271738 :   return gerepileuptoint(av,y);
    2493             : }
    2494             : 
    2495             : static GEN
    2496    14216238 : _Fp_mul(void *E, GEN x, GEN y) { return Fp_mul(x,y,(GEN)E); }
    2497             : static GEN
    2498     8134365 : _Fp_sqr(void *E, GEN x) { return Fp_sqr(x,(GEN)E); }
    2499             : static GEN
    2500       57564 : _Fp_one(void *E) { (void) E; return gen_1; }
    2501             : 
    2502             : GEN
    2503          77 : Fp_pow_init(GEN x, GEN n, long k, GEN p)
    2504          77 : { return gen_pow_init(x, n, k, (void*)p, &_Fp_sqr, &_Fp_mul); }
    2505             : 
    2506             : GEN
    2507       54096 : Fp_pow_table(GEN R, GEN n, GEN p)
    2508       54096 : { return gen_pow_table(R, n, (void*)p, &_Fp_one, &_Fp_mul); }
    2509             : 
    2510             : GEN
    2511        5924 : Fp_powers(GEN x, long n, GEN p)
    2512             : {
    2513        5924 :   if (lgefint(p) == 3)
    2514        2456 :     return Flv_to_ZV(Fl_powers(umodiu(x, uel(p, 2)), n, uel(p, 2)));
    2515        3468 :   return gen_powers(x, n, 1, (void*)p, _Fp_sqr, _Fp_mul, _Fp_one);
    2516             : }
    2517             : 
    2518             : GEN
    2519         497 : FpV_prod(GEN V, GEN p) { return gen_product(V, (void *)p, &_Fp_mul); }
    2520             : 
    2521             : static GEN
    2522    27823939 : _Fp_pow(void *E, GEN x, GEN n) { return Fp_pow(x,n,(GEN)E); }
    2523             : static GEN
    2524         168 : _Fp_rand(void *E) { return addiu(randomi(subiu((GEN)E,1)),1); }
    2525             : 
    2526             : static GEN Fp_easylog(void *E, GEN a, GEN g, GEN ord);
    2527             : static const struct bb_group Fp_star={_Fp_mul,_Fp_pow,_Fp_rand,hash_GEN,
    2528             :                                       equalii,equali1,Fp_easylog};
    2529             : 
    2530             : static GEN
    2531      927279 : _Fp_red(void *E, GEN x) { return Fp_red(x, (GEN)E); }
    2532             : static GEN
    2533     1221902 : _Fp_add(void *E, GEN x, GEN y) { (void) E; return addii(x,y); }
    2534             : static GEN
    2535     1123198 : _Fp_neg(void *E, GEN x) { (void) E; return negi(x); }
    2536             : static GEN
    2537      609346 : _Fp_rmul(void *E, GEN x, GEN y) { (void) E; return mulii(x,y); }
    2538             : static GEN
    2539       36641 : _Fp_inv(void *E, GEN x) { return Fp_inv(x,(GEN)E); }
    2540             : static int
    2541      265821 : _Fp_equal0(GEN x) { return signe(x)==0; }
    2542             : static GEN
    2543       33795 : _Fp_s(void *E, long x) { (void) E; return stoi(x); }
    2544             : 
    2545             : static const struct bb_field Fp_field={_Fp_red,_Fp_add,_Fp_rmul,_Fp_neg,
    2546             :                                         _Fp_inv,_Fp_equal0,_Fp_s};
    2547             : 
    2548        7868 : const struct bb_field *get_Fp_field(void **E, GEN p)
    2549        7868 : { *E = (void*)p; return &Fp_field; }
    2550             : 
    2551             : /*********************************************************************/
    2552             : /**               ORDER of INTEGERMOD x  in  (Z/nZ)*                **/
    2553             : /*********************************************************************/
    2554             : ulong
    2555      481766 : Fl_order(ulong a, ulong o, ulong p)
    2556             : {
    2557      481766 :   pari_sp av = avma;
    2558             :   GEN m, P, E;
    2559             :   long i;
    2560      481766 :   if (a==1) return 1;
    2561      438866 :   if (!o) o = p-1;
    2562      438866 :   m = factoru(o);
    2563      438866 :   P = gel(m,1);
    2564      438866 :   E = gel(m,2);
    2565     1251386 :   for (i = lg(P)-1; i; i--)
    2566             :   {
    2567      812520 :     ulong j, l = P[i], e = E[i], t = o / upowuu(l,e), y = Fl_powu(a, t, p);
    2568      812520 :     if (y == 1) o = t;
    2569      773077 :     else for (j = 1; j < e; j++)
    2570             :     {
    2571      385687 :       y = Fl_powu(y, l, p);
    2572      385687 :       if (y == 1) { o = t *  upowuu(l, j); break; }
    2573             :     }
    2574             :   }
    2575      438866 :   return gc_ulong(av, o);
    2576             : }
    2577             : 
    2578             : /*Find the exact order of a assuming a^o==1*/
    2579             : GEN
    2580       75380 : Fp_order(GEN a, GEN o, GEN p) {
    2581       75380 :   if (lgefint(p) == 3 && (!o || typ(o) == t_INT))
    2582             :   {
    2583          21 :     ulong pp = p[2], oo = (o && lgefint(o)==3)? uel(o,2): pp-1;
    2584          21 :     return utoi( Fl_order(umodiu(a, pp), oo, pp) );
    2585             :   }
    2586       75359 :   return gen_order(a, o, (void*)p, &Fp_star);
    2587             : }
    2588             : GEN
    2589          70 : Fp_factored_order(GEN a, GEN o, GEN p)
    2590          70 : { return gen_factored_order(a, o, (void*)p, &Fp_star); }
    2591             : 
    2592             : /* return order of a mod p^e, e > 0, pe = p^e */
    2593             : static GEN
    2594          70 : Zp_order(GEN a, GEN p, long e, GEN pe)
    2595             : {
    2596             :   GEN ap, op;
    2597          70 :   if (absequaliu(p, 2))
    2598             :   {
    2599          56 :     if (e == 1) return gen_1;
    2600          56 :     if (e == 2) return mod4(a) == 1? gen_1: gen_2;
    2601          49 :     if (mod4(a) == 1) op = gen_1; else { op = gen_2; a = Fp_sqr(a, pe); }
    2602             :   } else {
    2603          14 :     ap = (e == 1)? a: remii(a,p);
    2604          14 :     op = Fp_order(ap, subiu(p,1), p);
    2605          14 :     if (e == 1) return op;
    2606           0 :     a = Fp_pow(a, op, pe); /* 1 mod p */
    2607             :   }
    2608          49 :   if (equali1(a)) return op;
    2609           7 :   return mulii(op, powiu(p, e - Z_pval(subiu(a,1), p)));
    2610             : }
    2611             : 
    2612             : GEN
    2613          63 : znorder(GEN x, GEN o)
    2614             : {
    2615          63 :   pari_sp av = avma;
    2616             :   GEN b, a;
    2617             : 
    2618          63 :   if (typ(x) != t_INTMOD) pari_err_TYPE("znorder [t_INTMOD expected]",x);
    2619          56 :   b = gel(x,1); a = gel(x,2);
    2620          56 :   if (!equali1(gcdii(a,b))) pari_err_COPRIME("znorder", a,b);
    2621          49 :   if (!o)
    2622             :   {
    2623          35 :     GEN fa = Z_factor(b), P = gel(fa,1), E = gel(fa,2);
    2624          35 :     long i, l = lg(P);
    2625          35 :     o = gen_1;
    2626          70 :     for (i = 1; i < l; i++)
    2627             :     {
    2628          35 :       GEN p = gel(P,i);
    2629          35 :       long e = itos(gel(E,i));
    2630             : 
    2631          35 :       if (l == 2)
    2632          35 :         o = Zp_order(a, p, e, b);
    2633             :       else {
    2634           0 :         GEN pe = powiu(p,e);
    2635           0 :         o = lcmii(o, Zp_order(remii(a,pe), p, e, pe));
    2636             :       }
    2637             :     }
    2638          35 :     return gerepileuptoint(av, o);
    2639             :   }
    2640          14 :   return Fp_order(a, o, b);
    2641             : }
    2642             : 
    2643             : /*********************************************************************/
    2644             : /**               DISCRETE LOGARITHM  in  (Z/nZ)*                   **/
    2645             : /*********************************************************************/
    2646             : static GEN
    2647       58544 : Fp_log_halfgcd(ulong bnd, GEN C, GEN g, GEN p)
    2648             : {
    2649       58544 :   pari_sp av = avma;
    2650             :   GEN h1, h2, F, G;
    2651       58544 :   if (!Fp_ratlift(g,p,C,shifti(C,-1),&h1,&h2)) return gc_NULL(av);
    2652       35145 :   if ((F = Z_issmooth_fact(h1, bnd)) && (G = Z_issmooth_fact(h2, bnd)))
    2653             :   {
    2654         126 :     GEN M = cgetg(3, t_MAT);
    2655         126 :     gel(M,1) = vecsmall_concat(gel(F, 1),gel(G, 1));
    2656         126 :     gel(M,2) = vecsmall_concat(gel(F, 2),zv_neg_inplace(gel(G, 2)));
    2657         126 :     return gerepileupto(av, M);
    2658             :   }
    2659       35019 :   return gc_NULL(av);
    2660             : }
    2661             : 
    2662             : static GEN
    2663       58544 : Fp_log_find_rel(GEN b, ulong bnd, GEN C, GEN p, GEN *g, long *e)
    2664             : {
    2665             :   GEN rel;
    2666       58544 :   do { (*e)++; *g = Fp_mul(*g, b, p); rel = Fp_log_halfgcd(bnd, C, *g, p); }
    2667       58544 :   while (!rel);
    2668         126 :   return rel;
    2669             : }
    2670             : 
    2671             : struct Fp_log_rel
    2672             : {
    2673             :   GEN rel;
    2674             :   ulong prmax;
    2675             :   long nbrel, nbmax, nbgen;
    2676             : };
    2677             : 
    2678             : /* add u^e */
    2679             : static void
    2680        3157 : addifsmooth1(struct Fp_log_rel *r, GEN z, long u, long e)
    2681             : {
    2682        3157 :   pari_sp av = avma;
    2683        3157 :   long off = r->prmax+1;
    2684        3157 :   GEN F = cgetg(3, t_MAT);
    2685        3157 :   gel(F,1) = vecsmall_append(gel(z,1), off+u);
    2686        3157 :   gel(F,2) = vecsmall_append(gel(z,2), e);
    2687        3157 :   gel(r->rel,++r->nbrel) = gerepileupto(av, F);
    2688        3157 : }
    2689             : 
    2690             : /* add u^-1 v^-1 */
    2691             : static void
    2692      104083 : addifsmooth2(struct Fp_log_rel *r, GEN z, long u, long v)
    2693             : {
    2694      104083 :   pari_sp av = avma;
    2695      104083 :   long off = r->prmax+1;
    2696      104083 :   GEN P = mkvecsmall2(off+u,off+v), E = mkvecsmall2(-1,-1);
    2697      104083 :   GEN F = cgetg(3, t_MAT);
    2698      104083 :   gel(F,1) = vecsmall_concat(gel(z,1), P);
    2699      104083 :   gel(F,2) = vecsmall_concat(gel(z,2), E);
    2700      104083 :   gel(r->rel,++r->nbrel) = gerepileupto(av, F);
    2701      104083 : }
    2702             : 
    2703             : /*
    2704             : Let p=C^2+c
    2705             : Solve h = (C+x)*(C+a)-p = 0 [mod l]
    2706             : h= -c+x*(C+a)+C*a = 0  [mod l]
    2707             : x = (c-C*a)/(C+a) [mod l]
    2708             : h = -c+C*(x+a)+a*x
    2709             : */
    2710             : 
    2711             : GEN
    2712       40870 : Fp_log_sieve_worker(long a, long prmax, GEN C, GEN c, GEN Ci, GEN ci, GEN pi, GEN sz)
    2713             : {
    2714       40870 :   pari_sp ltop = avma;
    2715       40870 :   long i, j, th, n = lg(pi)-1, rel = 1;
    2716       40870 :   GEN sieve = zero_zv(a+2)+1;
    2717       40881 :   GEN L = cgetg(1+a+2, t_VEC);
    2718       40878 :   pari_sp av = avma;
    2719       40878 :   GEN z, h = addis(C,a);
    2720       40871 :   if ((z = Z_issmooth_fact(h, prmax)))
    2721             :   {
    2722        3007 :     gel(L, rel++) = mkvec2(z, mkvecsmall3(1, a, -1));
    2723        3012 :     av = avma;
    2724             :   }
    2725    16733507 :   for (i=1; i<=n; i++)
    2726             :   {
    2727    16691568 :     ulong li = pi[i], s = sz[i], al = a % li;
    2728    16691568 :     ulong u, iv = Fl_invsafe(Fl_add(Ci[i],al,li),li);
    2729    16940464 :     if (!iv) continue;
    2730    16516293 :     u = Fl_mul(Fl_sub(ci[i],Fl_mul(Ci[i],al,li),li), iv ,li);
    2731    77693542 :     for(j = u; j<=a; j+=li) sieve[j] += s;
    2732             :   }
    2733       41939 :   if (a)
    2734             :   {
    2735       40738 :     long e = expi(mulis(C,a));
    2736       40786 :     th = e - expu(e) - 1;
    2737        1201 :   } else th = -1;
    2738    28019883 :   for (j=0; j<a; j++)
    2739    27978475 :     if (sieve[j]>=th)
    2740             :     {
    2741      119404 :       GEN h = addiu(subii(muliu(C,a+j),c), a*j);
    2742      119434 :       if ((z = Z_issmooth_fact(h, prmax)))
    2743             :       {
    2744      109720 :         gel(L, rel++) = mkvec2(z, mkvecsmall3(2, a, j));
    2745      109927 :         av = avma;
    2746        9409 :       } else set_avma(av);
    2747             :     }
    2748             :   /* j = a */
    2749       41408 :   if (sieve[a]>=th)
    2750             :   {
    2751         476 :     GEN h = addiu(subii(muliu(C,2*a),c), a*a);
    2752         476 :     if ((z = Z_issmooth_fact(h, prmax)))
    2753         385 :       gel(L, rel++) = mkvec2(z, mkvecsmall3(1, a, -2));
    2754             :   }
    2755       41408 :   setlg(L, rel); return gerepilecopy(ltop, L);
    2756             : }
    2757             : 
    2758             : static long
    2759          63 : Fp_log_sieve(struct Fp_log_rel *r, GEN C, GEN c, GEN Ci, GEN ci, GEN pi, GEN sz)
    2760             : {
    2761             :   struct pari_mt pt;
    2762          63 :   long i, j, nb = 0;
    2763          63 :   GEN worker = snm_closure(is_entry("_Fp_log_sieve_worker"),
    2764             :                mkvecn(7, utoi(r->prmax), C, c, Ci, ci, pi, sz));
    2765          63 :   long running, pending = 0;
    2766          63 :   GEN W = zerovec(r->nbgen);
    2767          63 :   mt_queue_start_lim(&pt, worker, r->nbgen);
    2768       41085 :   for (i = 0; (running = (i < r->nbgen)) || pending; i++)
    2769             :   {
    2770             :     GEN done;
    2771             :     long idx;
    2772       41022 :     mt_queue_submit(&pt, i, running ? mkvec(stoi(i)): NULL);
    2773       41022 :     done = mt_queue_get(&pt, &idx, &pending);
    2774       41022 :     if (!done || lg(done)==1) continue;
    2775       35917 :     gel(W, idx+1) = done;
    2776       35917 :     nb += lg(done)-1;
    2777       35917 :     if (DEBUGLEVEL && (i&127)==0)
    2778           0 :       err_printf("%ld%% ",100*nb/r->nbmax);
    2779             :   }
    2780          63 :   mt_queue_end(&pt);
    2781       39550 :   for(j = 1; j <= r->nbgen && r->nbrel < r->nbmax; j++)
    2782             :   {
    2783             :     long ll, m;
    2784       39487 :     GEN L = gel(W,j);
    2785       39487 :     if (isintzero(L)) continue;
    2786       34531 :     ll = lg(L);
    2787      141771 :     for (m=1; m<ll && r->nbrel < r->nbmax ; m++)
    2788             :     {
    2789      107240 :       GEN Lm = gel(L,m), h = gel(Lm, 1), v = gel(Lm, 2);
    2790      107240 :       if (v[1] == 1)
    2791        3157 :         addifsmooth1(r, h, v[2], v[3]);
    2792             :       else
    2793      104083 :         addifsmooth2(r, h, v[2], v[3]);
    2794             :     }
    2795             :   }
    2796          63 :   return j;
    2797             : }
    2798             : 
    2799             : static GEN
    2800         735 : ECP_psi(GEN x, GEN y)
    2801             : {
    2802         735 :   long prec = realprec(x);
    2803         735 :   GEN lx = glog(x, prec), ly = glog(y, prec);
    2804         735 :   GEN u = gdiv(lx, ly);
    2805         735 :   return gpow(u, gneg(u),prec);
    2806             : }
    2807             : 
    2808             : struct computeG
    2809             : {
    2810             :   GEN C;
    2811             :   long bnd, nbi;
    2812             : };
    2813             : 
    2814             : static GEN
    2815         735 : _computeG(void *E, GEN gen)
    2816             : {
    2817         735 :   struct computeG * d = (struct computeG *) E;
    2818         735 :   GEN ps = ECP_psi(gmul(gen,d->C), stoi(d->bnd));
    2819         735 :   return gsub(gmul(gsqr(gen),ps),gmul2n(gaddgs(gen,d->nbi),2));
    2820             : }
    2821             : 
    2822             : static long
    2823          63 : compute_nbgen(GEN C, long bnd, long nbi)
    2824             : {
    2825             :   struct computeG d;
    2826          63 :   d.C = shifti(C, 1);
    2827          63 :   d.bnd = bnd;
    2828          63 :   d.nbi = nbi;
    2829          63 :   return itos(ground(zbrent((void*)&d, _computeG, gen_2, stoi(bnd), DEFAULTPREC)));
    2830             : }
    2831             : 
    2832             : static GEN
    2833        1714 : _psi(void*E, GEN y)
    2834             : {
    2835        1714 :   GEN lx = (GEN) E;
    2836        1714 :   long prec = realprec(lx);
    2837        1714 :   GEN ly = glog(y, prec);
    2838        1714 :   GEN u = gdiv(lx, ly);
    2839        1714 :   return gsub(gdiv(y ,ly), gpow(u, u, prec));
    2840             : }
    2841             : 
    2842             : static GEN
    2843          63 : opt_param(GEN x, long prec)
    2844             : {
    2845          63 :   return zbrent((void*)glog(x,prec), _psi, gen_2, x, prec);
    2846             : }
    2847             : 
    2848             : static GEN
    2849          63 : check_kernel(long nbg, long N, long prmax, GEN C, GEN M, GEN p, GEN m)
    2850             : {
    2851          63 :   pari_sp av = avma;
    2852          63 :   long lM = lg(M)-1, nbcol = lM;
    2853          63 :   long tbs = maxss(1, expu(nbg/expi(m)));
    2854             :   for (;;)
    2855          14 :   {
    2856          77 :     GEN K = FpMs_leftkernel_elt_col(M, nbcol, N, m);
    2857             :     GEN tab;
    2858          77 :     long i, f=0;
    2859          77 :     long l = lg(K), lm = lgefint(m);
    2860          77 :     GEN idx = diviiexact(subiu(p,1),m), g;
    2861             :     pari_timer ti;
    2862          77 :     if (DEBUGLEVEL) timer_start(&ti);
    2863         154 :     for(i=1; i<l; i++)
    2864         154 :       if (signe(gel(K,i)))
    2865          77 :         break;
    2866          77 :     g = Fp_pow(utoi(i), idx, p);
    2867          77 :     tab = Fp_pow_init(g, p, tbs, p);
    2868          77 :     K = FpC_Fp_mul(K, Fp_inv(gel(K,i), m), m);
    2869      128464 :     for(i=1; i<l; i++)
    2870             :     {
    2871      128387 :       GEN k = gel(K,i);
    2872      128387 :       GEN j = i<=prmax ? utoi(i): addis(C,i-(prmax+1));
    2873      128387 :       if (signe(k)==0 || !equalii(Fp_pow_table(tab, k, p), Fp_pow(j, idx, p)))
    2874       76391 :         gel(K,i) = cgetineg(lm);
    2875             :       else
    2876       51996 :         f++;
    2877             :     }
    2878          77 :     if (DEBUGLEVEL) timer_printf(&ti,"found %ld/%ld logs", f, nbg);
    2879          77 :     if(f > (nbg>>1)) return gerepileupto(av, K);
    2880        4585 :     for(i=1; i<=nbcol; i++)
    2881             :     {
    2882        4571 :       long a = 1+random_Fl(lM);
    2883        4571 :       swap(gel(M,a),gel(M,i));
    2884             :     }
    2885          14 :     if (4*nbcol>5*nbg) nbcol = nbcol*9/10;
    2886             :   }
    2887             : }
    2888             : 
    2889             : static GEN
    2890         126 : Fp_log_find_ind(GEN a, GEN K, long prmax, GEN C, GEN p, GEN m)
    2891             : {
    2892         126 :   pari_sp av=avma;
    2893         126 :   GEN aa = gen_1;
    2894         126 :   long AV = 0;
    2895             :   for(;;)
    2896           0 :   {
    2897         126 :     GEN A = Fp_log_find_rel(a, prmax, C, p, &aa, &AV);
    2898         126 :     GEN F = gel(A,1), E = gel(A,2);
    2899         126 :     GEN Ao = gen_0;
    2900         126 :     long i, l = lg(F);
    2901         811 :     for(i=1; i<l; i++)
    2902             :     {
    2903         685 :       GEN Ki = gel(K,F[i]);
    2904         685 :       if (signe(Ki)<0) break;
    2905         685 :       Ao = addii(Ao, mulis(Ki, E[i]));
    2906             :     }
    2907         126 :     if (i==l) return Fp_divu(Ao, AV, m);
    2908           0 :     aa = gerepileuptoint(av, aa);
    2909             :   }
    2910             : }
    2911             : 
    2912             : static GEN
    2913          63 : Fp_log_index(GEN a, GEN b, GEN m, GEN p)
    2914             : {
    2915          63 :   pari_sp av = avma, av2;
    2916          63 :   long i, j, nbi, nbr = 0, nbrow, nbg;
    2917             :   GEN C, c, Ci, ci, pi, pr, sz, l, Ao, Bo, K, d, p_1;
    2918             :   pari_timer ti;
    2919             :   struct Fp_log_rel r;
    2920          63 :   ulong bnds = itou(roundr_safe(opt_param(sqrti(p),DEFAULTPREC)));
    2921          63 :   ulong bnd = 4*bnds;
    2922          63 :   if (!bnds || cmpii(sqru(bnds),m)>=0) return NULL;
    2923             : 
    2924          63 :   p_1 = subiu(p,1);
    2925          63 :   if (!is_pm1(gcdii(m,diviiexact(p_1,m))))
    2926           0 :     m = diviiexact(p_1, Z_ppo(p_1, m));
    2927          63 :   pr = primes_upto_zv(bnd);
    2928          63 :   nbi = lg(pr)-1;
    2929          63 :   C = sqrtremi(p, &c);
    2930          63 :   av2 = avma;
    2931       12796 :   for (i = 1; i <= nbi; ++i)
    2932             :   {
    2933       12733 :     ulong lp = pr[i];
    2934       26894 :     while (lp <= bnd)
    2935             :     {
    2936       14161 :       nbr++;
    2937       14161 :       lp *= pr[i];
    2938             :     }
    2939             :   }
    2940          63 :   pi = cgetg(nbr+1,t_VECSMALL);
    2941          63 :   Ci = cgetg(nbr+1,t_VECSMALL);
    2942          63 :   ci = cgetg(nbr+1,t_VECSMALL);
    2943          63 :   sz = cgetg(nbr+1,t_VECSMALL);
    2944       12796 :   for (i = 1, j = 1; i <= nbi; ++i)
    2945             :   {
    2946       12733 :     ulong lp = pr[i], sp = expu(2*lp-1);
    2947       26894 :     while (lp <= bnd)
    2948             :     {
    2949       14161 :       pi[j] = lp;
    2950       14161 :       Ci[j] = umodiu(C, lp);
    2951       14161 :       ci[j] = umodiu(c, lp);
    2952       14161 :       sz[j] = sp;
    2953       14161 :       lp *= pr[i];
    2954       14161 :       j++;
    2955             :     }
    2956             :   }
    2957          63 :   r.nbrel = 0;
    2958          63 :   r.nbgen = compute_nbgen(C, bnd, nbi);
    2959          63 :   r.nbmax = 2*(nbi+r.nbgen);
    2960          63 :   r.rel = cgetg(r.nbmax+1,t_VEC);
    2961          63 :   r.prmax = pr[nbi];
    2962          63 :   if (DEBUGLEVEL)
    2963             :   {
    2964           0 :     err_printf("bnd=%lu Size FB=%ld extra gen=%ld \n", bnd, nbi, r.nbgen);
    2965           0 :     timer_start(&ti);
    2966             :   }
    2967          63 :   nbg = Fp_log_sieve(&r, C, c, Ci, ci, pi, sz);
    2968          63 :   nbrow = r.prmax + nbg;
    2969          63 :   if (DEBUGLEVEL)
    2970             :   {
    2971           0 :     err_printf("\n");
    2972           0 :     timer_printf(&ti," %ld relations, %ld generators", r.nbrel, nbi+nbg);
    2973             :   }
    2974          63 :   setlg(r.rel,r.nbrel+1);
    2975          63 :   r.rel = gerepilecopy(av2, r.rel);
    2976          63 :   K = check_kernel(nbi+nbrow-r.prmax, nbrow, r.prmax, C, r.rel, p, m);
    2977          63 :   if (DEBUGLEVEL) timer_start(&ti);
    2978          63 :   Ao = Fp_log_find_ind(a, K, r.prmax, C, p, m);
    2979          63 :   if (DEBUGLEVEL) timer_printf(&ti," log element");
    2980          63 :   Bo = Fp_log_find_ind(b, K, r.prmax, C, p, m);
    2981          63 :   if (DEBUGLEVEL) timer_printf(&ti," log generator");
    2982          63 :   d = gcdii(Ao,Bo);
    2983          63 :   l = Fp_div(diviiexact(Ao, d) ,diviiexact(Bo, d), m);
    2984          63 :   if (!equalii(a,Fp_pow(b,l,p))) pari_err_BUG("Fp_log_index");
    2985          63 :   return gerepileuptoint(av, l);
    2986             : }
    2987             : 
    2988             : static int
    2989     4652477 : Fp_log_use_index(long e, long p)
    2990             : {
    2991     4652477 :   return (e >= 27 && 20*(p+6)<=e*e);
    2992             : }
    2993             : 
    2994             : /* Trivial cases a = 1, -1. Return x s.t. g^x = a or [] if no such x exist */
    2995             : static GEN
    2996     8447249 : Fp_easylog(void *E, GEN a, GEN g, GEN ord)
    2997             : {
    2998     8447249 :   pari_sp av = avma;
    2999     8447249 :   GEN p = (GEN)E;
    3000             :   /* assume a reduced mod p, p not necessarily prime */
    3001     8447249 :   if (equali1(a)) return gen_0;
    3002             :   /* p > 2 */
    3003     5427404 :   if (equalii(subiu(p,1), a))  /* -1 */
    3004             :   {
    3005             :     pari_sp av2;
    3006             :     GEN t;
    3007     1318101 :     ord = get_arith_Z(ord);
    3008     1318101 :     if (mpodd(ord)) { set_avma(av); return cgetg(1, t_VEC); } /* no solution */
    3009     1318087 :     t = shifti(ord,-1); /* only possible solution */
    3010     1318087 :     av2 = avma;
    3011     1318087 :     if (!equalii(Fp_pow(g, t, p), a)) { set_avma(av); return cgetg(1, t_VEC); }
    3012     1318059 :     set_avma(av2); return gerepileuptoint(av, t);
    3013             :   }
    3014     4109308 :   if (typ(ord)==t_INT && BPSW_psp(p) && Fp_log_use_index(expi(ord),expi(p)))
    3015          63 :     return Fp_log_index(a, g, ord, p);
    3016     4109245 :   return gc_NULL(av); /* not easy */
    3017             : }
    3018             : 
    3019             : GEN
    3020     3917819 : Fp_log(GEN a, GEN g, GEN ord, GEN p)
    3021             : {
    3022     3917819 :   GEN v = get_arith_ZZM(ord);
    3023     3917798 :   GEN F = gmael(v,2,1);
    3024     3917798 :   long lF = lg(F)-1, lmax;
    3025     3917798 :   if (lF == 0) return equali1(a)? gen_0: cgetg(1, t_VEC);
    3026     3917770 :   lmax = expi(gel(F,lF));
    3027     3917768 :   if (BPSW_psp(p) && Fp_log_use_index(lmax,expi(p)))
    3028          91 :     v = mkvec2(gel(v,1),ZM_famat_limit(gel(v,2),int2n(27)));
    3029     3917768 :   return gen_PH_log(a,g,v,(void*)p,&Fp_star);
    3030             : }
    3031             : 
    3032             : /* assume !(p & HIGHMASK) */
    3033             : static ulong
    3034      131786 : Fl_log_naive(ulong a, ulong g, ulong ord, ulong p)
    3035             : {
    3036      131786 :   ulong i, h=1;
    3037      361301 :   for (i = 0; i < ord; i++, h = (h * g) % p)
    3038      361302 :     if (a==h) return i;
    3039           0 :   return ~0UL;
    3040             : }
    3041             : 
    3042             : static ulong
    3043       24558 : Fl_log_naive_pre(ulong a, ulong g, ulong ord, ulong p, ulong pi)
    3044             : {
    3045       24558 :   ulong i, h=1;
    3046       62902 :   for (i = 0; i < ord; i++, h = Fl_mul_pre(h, g, p, pi))
    3047       62902 :     if (a==h) return i;
    3048           0 :   return ~0UL;
    3049             : }
    3050             : 
    3051             : static ulong
    3052           0 : Fl_log_Fp(ulong a, ulong g, ulong ord, ulong p)
    3053             : {
    3054           0 :   pari_sp av = avma;
    3055           0 :   GEN r = Fp_log(utoi(a),utoi(g),utoi(ord),utoi(p));
    3056           0 :   return gc_ulong(av, typ(r)==t_INT ? itou(r): ~0UL);
    3057             : }
    3058             : 
    3059             : /* allow pi = 0 */
    3060             : ulong
    3061       24957 : Fl_log_pre(ulong a, ulong g, ulong ord, ulong p, ulong pi)
    3062             : {
    3063       24957 :   if (!pi) return Fl_log(a, g, ord, p);
    3064       24558 :   if (ord <= 200) return Fl_log_naive_pre(a, g, ord, p, pi);
    3065           0 :   return Fl_log_Fp(a, g, ord, p);
    3066             : }
    3067             : 
    3068             : ulong
    3069      131786 : Fl_log(ulong a, ulong g, ulong ord, ulong p)
    3070             : {
    3071      131786 :   if (ord <= 200)
    3072           0 :     return (p&HIGHMASK)? Fl_log_naive_pre(a, g, ord, p, get_Fl_red(p))
    3073      131786 :                        : Fl_log_naive(a, g, ord, p);
    3074           0 :   return Fl_log_Fp(a, g, ord, p);
    3075             : }
    3076             : 
    3077             : /* find x such that h = g^x mod N > 1, N = prod_{i <= l} P[i]^E[i], P[i] prime.
    3078             :  * PHI[l] = eulerphi(N / P[l]^E[l]).   Destroys P/E */
    3079             : static GEN
    3080         126 : znlog_rec(GEN h, GEN g, GEN N, GEN P, GEN E, GEN PHI)
    3081             : {
    3082         126 :   long l = lg(P) - 1, e = E[l];
    3083         126 :   GEN p = gel(P, l), phi = gel(PHI,l), pe = e == 1? p: powiu(p, e);
    3084             :   GEN a,b, hp,gp, hpe,gpe, ogpe; /* = order(g mod p^e) | p^(e-1)(p-1) */
    3085             : 
    3086         126 :   if (l == 1) {
    3087          98 :     hpe = h;
    3088          98 :     gpe = g;
    3089             :   } else {
    3090          28 :     hpe = modii(h, pe);
    3091          28 :     gpe = modii(g, pe);
    3092             :   }
    3093         126 :   if (e == 1) {
    3094          42 :     hp = hpe;
    3095          42 :     gp = gpe;
    3096             :   } else {
    3097          84 :     hp = remii(hpe, p);
    3098          84 :     gp = remii(gpe, p);
    3099             :   }
    3100         126 :   if (hp == gen_0 || gp == gen_0) return NULL;
    3101         105 :   if (absequaliu(p, 2))
    3102             :   {
    3103          35 :     GEN n = int2n(e);
    3104          35 :     ogpe = Zp_order(gpe, gen_2, e, n);
    3105          35 :     a = Fp_log(hpe, gpe, ogpe, n);
    3106          35 :     if (typ(a) != t_INT) return NULL;
    3107             :   }
    3108             :   else
    3109             :   { /* Avoid black box groups: (Z/p^2)^* / (Z/p)^* ~ (Z/pZ, +), where DL
    3110             :        is trivial */
    3111             :     /* [order(gp), factor(order(gp))] */
    3112          70 :     GEN v = Fp_factored_order(gp, subiu(p,1), p);
    3113          70 :     GEN ogp = gel(v,1);
    3114          70 :     if (!equali1(Fp_pow(hp, ogp, p))) return NULL;
    3115          70 :     a = Fp_log(hp, gp, v, p);
    3116          70 :     if (typ(a) != t_INT) return NULL;
    3117          70 :     if (e == 1) ogpe = ogp;
    3118             :     else
    3119             :     { /* find a s.t. g^a = h (mod p^e), p odd prime, e > 0, (h,p) = 1 */
    3120             :       /* use p-adic log: O(log p + e) mul*/
    3121             :       long vpogpe, vpohpe;
    3122             : 
    3123          28 :       hpe = Fp_mul(hpe, Fp_pow(gpe, negi(a), pe), pe);
    3124          28 :       gpe = Fp_pow(gpe, ogp, pe);
    3125             :       /* g,h = 1 mod p; compute b s.t. h = g^b */
    3126             : 
    3127             :       /* v_p(order g mod pe) */
    3128          28 :       vpogpe = equali1(gpe)? 0: e - Z_pval(subiu(gpe,1), p);
    3129             :       /* v_p(order h mod pe) */
    3130          28 :       vpohpe = equali1(hpe)? 0: e - Z_pval(subiu(hpe,1), p);
    3131          28 :       if (vpohpe > vpogpe) return NULL;
    3132             : 
    3133          28 :       ogpe = mulii(ogp, powiu(p, vpogpe)); /* order g mod p^e */
    3134          28 :       if (is_pm1(gpe)) return is_pm1(hpe)? a: NULL;
    3135          28 :       b = gdiv(Qp_log(cvtop(hpe, p, e)), Qp_log(cvtop(gpe, p, e)));
    3136          28 :       a = addii(a, mulii(ogp, padic_to_Q(b)));
    3137             :     }
    3138             :   }
    3139             :   /* gp^a = hp => x = a mod ogpe => generalized Pohlig-Hellman strategy */
    3140          91 :   if (l == 1) return a;
    3141             : 
    3142          28 :   N = diviiexact(N, pe); /* make N coprime to p */
    3143          28 :   h = Fp_mul(h, Fp_pow(g, modii(negi(a), phi), N), N);
    3144          28 :   g = Fp_pow(g, modii(ogpe, phi), N);
    3145          28 :   setlg(P, l); /* remove last element */
    3146          28 :   setlg(E, l);
    3147          28 :   b = znlog_rec(h, g, N, P, E, PHI);
    3148          28 :   if (!b) return NULL;
    3149          28 :   return addmulii(a, b, ogpe);
    3150             : }
    3151             : 
    3152             : static GEN
    3153          98 : get_PHI(GEN P, GEN E)
    3154             : {
    3155          98 :   long i, l = lg(P);
    3156          98 :   GEN PHI = cgetg(l, t_VEC);
    3157          98 :   gel(PHI,1) = gen_1;
    3158         126 :   for (i=1; i<l-1; i++)
    3159             :   {
    3160          28 :     GEN t, p = gel(P,i);
    3161          28 :     long e = E[i];
    3162          28 :     t = mulii(powiu(p, e-1), subiu(p,1));
    3163          28 :     if (i > 1) t = mulii(t, gel(PHI,i));
    3164          28 :     gel(PHI,i+1) = t;
    3165             :   }
    3166          98 :   return PHI;
    3167             : }
    3168             : 
    3169             : GEN
    3170         238 : znlog(GEN h, GEN g, GEN o)
    3171             : {
    3172         238 :   pari_sp av = avma;
    3173             :   GEN N, fa, P, E, x;
    3174         238 :   switch (typ(g))
    3175             :   {
    3176          28 :     case t_PADIC:
    3177             :     {
    3178          28 :       GEN p = gel(g,2);
    3179          28 :       long v = valp(g);
    3180          28 :       if (v < 0) pari_err_DIM("znlog");
    3181          28 :       if (v > 0) {
    3182           0 :         long k = gvaluation(h, p);
    3183           0 :         if (k % v) return cgetg(1,t_VEC);
    3184           0 :         k /= v;
    3185           0 :         if (!gequal(h, gpowgs(g,k))) { set_avma(av); return cgetg(1,t_VEC); }
    3186           0 :         return gc_stoi(av, k);
    3187             :       }
    3188          28 :       N = gel(g,3);
    3189          28 :       g = Rg_to_Fp(g, N);
    3190          28 :       break;
    3191             :     }
    3192         203 :     case t_INTMOD:
    3193         203 :       N = gel(g,1);
    3194         203 :       g = gel(g,2); break;
    3195           7 :     default: pari_err_TYPE("znlog", g);
    3196             :       return NULL; /* LCOV_EXCL_LINE */
    3197             :   }
    3198         231 :   if (equali1(N)) { set_avma(av); return gen_0; }
    3199         231 :   h = Rg_to_Fp(h, N);
    3200         224 :   if (o) return gerepileupto(av, Fp_log(h, g, o, N));
    3201          98 :   fa = Z_factor(N);
    3202          98 :   P = gel(fa,1);
    3203          98 :   E = vec_to_vecsmall(gel(fa,2));
    3204          98 :   x = znlog_rec(h, g, N, P, E, get_PHI(P,E));
    3205          98 :   if (!x) { set_avma(av); return cgetg(1,t_VEC); }
    3206          63 :   return gerepileuptoint(av, x);
    3207             : }
    3208             : 
    3209             : GEN
    3210      173306 : Fp_sqrtn(GEN a, GEN n, GEN p, GEN *zeta)
    3211             : {
    3212      173306 :   if (lgefint(p)==3)
    3213             :   {
    3214      172872 :     long nn = itos_or_0(n);
    3215      172872 :     if (nn)
    3216             :     {
    3217      172872 :       ulong pp = p[2];
    3218             :       ulong uz;
    3219      172872 :       ulong r = Fl_sqrtn(umodiu(a,pp),nn,pp, zeta ? &uz:NULL);
    3220      172851 :       if (r==ULONG_MAX) return NULL;
    3221      172795 :       if (zeta) *zeta = utoi(uz);
    3222      172795 :       return utoi(r);
    3223             :     }
    3224             :   }
    3225         434 :   a = modii(a,p);
    3226         434 :   if (!signe(a))
    3227             :   {
    3228           0 :     if (zeta) *zeta = gen_1;
    3229           0 :     if (signe(n) < 0) pari_err_INV("Fp_sqrtn", mkintmod(gen_0,p));
    3230           0 :     return gen_0;
    3231             :   }
    3232         434 :   if (absequaliu(n,2))
    3233             :   {
    3234         238 :     if (zeta) *zeta = subiu(p,1);
    3235         238 :     return signe(n) > 0 ? Fp_sqrt(a,p): Fp_sqrt(Fp_inv(a, p),p);
    3236             :   }
    3237         196 :   return gen_Shanks_sqrtn(a,n,subiu(p,1),zeta,(void*)p,&Fp_star);
    3238             : }
    3239             : 
    3240             : /*********************************************************************/
    3241             : /**                              FACTORIAL                          **/
    3242             : /*********************************************************************/
    3243             : GEN
    3244       90239 : mulu_interval_step(ulong a, ulong b, ulong step)
    3245             : {
    3246       90239 :   pari_sp av = avma;
    3247             :   ulong k, l, N, n;
    3248             :   long lx;
    3249             :   GEN x;
    3250             : 
    3251       90239 :   if (!a) return gen_0;
    3252       90239 :   if (step == 1) return mulu_interval(a, b);
    3253       90239 :   n = 1 + (b-a) / step;
    3254       90239 :   b -= (b-a) % step;
    3255       90239 :   if (n < 61)
    3256             :   {
    3257       88802 :     if (n == 1) return utoipos(a);
    3258       68316 :     x = muluu(a,a+step); if (n == 2) return x;
    3259      535893 :     for (k=a+2*step; k<=b; k+=step) x = mului(k,x);
    3260       53411 :     return gerepileuptoint(av, x);
    3261             :   }
    3262             :   /* step | b-a */
    3263        1437 :   lx = 1; x = cgetg(2 + n/2, t_VEC);
    3264        1440 :   N = b + a;
    3265        1440 :   for (k = a;; k += step)
    3266             :   {
    3267      229975 :     l = N - k; if (l <= k) break;
    3268      228535 :     gel(x,lx++) = muluu(k,l);
    3269             :   }
    3270        1440 :   if (l == k) gel(x,lx++) = utoipos(k);
    3271        1440 :   setlg(x, lx);
    3272        1440 :   return gerepileuptoint(av, ZV_prod(x));
    3273             : }
    3274             : /* return a * (a+1) * ... * b. Assume a <= b  [ note: factoring out powers of 2
    3275             :  * first is slower ... ] */
    3276             : GEN
    3277      158209 : mulu_interval(ulong a, ulong b)
    3278             : {
    3279      158209 :   pari_sp av = avma;
    3280             :   ulong k, l, N, n;
    3281             :   long lx;
    3282             :   GEN x;
    3283             : 
    3284      158209 :   if (!a) return gen_0;
    3285      158209 :   n = b - a + 1;
    3286      158209 :   if (n < 61)
    3287             :   {
    3288      157481 :     if (n == 1) return utoipos(a);
    3289      107165 :     x = muluu(a,a+1); if (n == 2) return x;
    3290      401123 :     for (k=a+2; k<b; k++) x = mului(k,x);
    3291             :     /* avoid k <= b: broken if b = ULONG_MAX */
    3292       93113 :     return gerepileuptoint(av, mului(b,x));
    3293             :   }
    3294         728 :   lx = 1; x = cgetg(2 + n/2, t_VEC);
    3295         733 :   N = b + a;
    3296         733 :   for (k = a;; k++)
    3297             :   {
    3298       27226 :     l = N - k; if (l <= k) break;
    3299       26497 :     gel(x,lx++) = muluu(k,l);
    3300             :   }
    3301         729 :   if (l == k) gel(x,lx++) = utoipos(k);
    3302         729 :   setlg(x, lx);
    3303         730 :   return gerepileuptoint(av, ZV_prod(x));
    3304             : }
    3305             : GEN
    3306         560 : muls_interval(long a, long b)
    3307             : {
    3308         560 :   pari_sp av = avma;
    3309         560 :   long lx, k, l, N, n = b - a + 1;
    3310             :   GEN x;
    3311             : 
    3312         560 :   if (a <= 0 && b >= 0) return gen_0;
    3313         287 :   if (n < 61)
    3314             :   {
    3315         287 :     x = stoi(a);
    3316         511 :     for (k=a+1; k<=b; k++) x = mulsi(k,x);
    3317         287 :     return gerepileuptoint(av, x);
    3318             :   }
    3319           0 :   lx = 1; x = cgetg(2 + n/2, t_VEC);
    3320           0 :   N = b + a;
    3321           0 :   for (k = a;; k++)
    3322             :   {
    3323           0 :     l = N - k; if (l <= k) break;
    3324           0 :     gel(x,lx++) = mulss(k,l);
    3325             :   }
    3326           0 :   if (l == k) gel(x,lx++) = stoi(k);
    3327           0 :   setlg(x, lx);
    3328           0 :   return gerepileuptoint(av, ZV_prod(x));
    3329             : }
    3330             : 
    3331             : GEN
    3332         105 : mpprimorial(long n)
    3333             : {
    3334         105 :   pari_sp av = avma;
    3335         105 :   if (n <= 12) switch(n)
    3336             :   {
    3337          14 :     case 0: case 1: return gen_1;
    3338           7 :     case 2: return gen_2;
    3339          14 :     case 3: case 4: return utoipos(6);
    3340          14 :     case 5: case 6: return utoipos(30);
    3341          28 :     case 7: case 8: case 9: case 10: return utoipos(210);
    3342          14 :     case 11: case 12: return utoipos(2310);
    3343           7 :     default: pari_err_DOMAIN("primorial", "argument","<",gen_0,stoi(n));
    3344             :   }
    3345           7 :   return gerepileuptoint(av, zv_prod_Z(primes_upto_zv(n)));
    3346             : }
    3347             : 
    3348             : GEN
    3349      504942 : mpfact(long n)
    3350             : {
    3351      504942 :   pari_sp av = avma;
    3352             :   GEN a, v;
    3353             :   long k;
    3354      504942 :   if (n <= 12) switch(n)
    3355             :   {
    3356      424243 :     case 0: case 1: return gen_1;
    3357       37091 :     case 2: return gen_2;
    3358        3409 :     case 3: return utoipos(6);
    3359        4222 :     case 4: return utoipos(24);
    3360        2873 :     case 5: return utoipos(120);
    3361        2528 :     case 6: return utoipos(720);
    3362        2434 :     case 7: return utoipos(5040);
    3363        2423 :     case 8: return utoipos(40320);
    3364        2444 :     case 9: return utoipos(362880);
    3365        2708 :     case 10:return utoipos(3628800);
    3366        1437 :     case 11:return utoipos(39916800);
    3367         591 :     case 12:return utoipos(479001600);
    3368           0 :     default: pari_err_DOMAIN("factorial", "argument","<",gen_0,stoi(n));
    3369             :   }
    3370       18541 :   v = cgetg(expu(n) + 2, t_VEC);
    3371       18548 :   for (k = 1;; k++)
    3372       86421 :   {
    3373      104969 :     long m = n >> (k-1), l;
    3374      104969 :     if (m <= 2) break;
    3375       86428 :     l = (1 + (n >> k)) | 1;
    3376             :     /* product of odd numbers in ]n / 2^k, 2 / 2^(k-1)] */
    3377       86428 :     a = mulu_interval_step(l, m, 2);
    3378       86304 :     gel(v,k) = k == 1? a: powiu(a, k);
    3379             :   }
    3380       86434 :   a = gel(v,--k); while (--k) a = mulii(a, gel(v,k));
    3381       18534 :   a = shifti(a, factorial_lval(n, 2));
    3382       18533 :   return gerepileuptoint(av, a);
    3383             : }
    3384             : 
    3385             : ulong
    3386       56693 : factorial_Fl(long n, ulong p)
    3387             : {
    3388             :   long k;
    3389             :   ulong v;
    3390       56693 :   if (p <= (ulong)n) return 0;
    3391       56693 :   v = Fl_powu(2, factorial_lval(n, 2), p);
    3392       56767 :   for (k = 1;; k++)
    3393      142333 :   {
    3394      199100 :     long m = n >> (k-1), l, i;
    3395      199100 :     ulong a = 1;
    3396      199100 :     if (m <= 2) break;
    3397      142347 :     l = (1 + (n >> k)) | 1;
    3398             :     /* product of odd numbers in ]n / 2^k, 2 / 2^(k-1)] */
    3399      778667 :     for (i=l; i<=m; i+=2)
    3400      636323 :       a = Fl_mul(a, i, p);
    3401      142344 :     v = Fl_mul(v, k == 1? a: Fl_powu(a, k, p), p);
    3402             :   }
    3403       56753 :   return v;
    3404             : }
    3405             : 
    3406             : GEN
    3407         158 : factorial_Fp(long n, GEN p)
    3408             : {
    3409         158 :   pari_sp av = avma;
    3410             :   long k;
    3411         158 :   GEN v = Fp_powu(gen_2, factorial_lval(n, 2), p);
    3412         158 :   for (k = 1;; k++)
    3413         344 :   {
    3414         502 :     long m = n >> (k-1), l, i;
    3415         502 :     GEN a = gen_1;
    3416         502 :     if (m <= 2) break;
    3417         344 :     l = (1 + (n >> k)) | 1;
    3418             :     /* product of odd numbers in ]n / 2^k, 2 / 2^(k-1)] */
    3419         850 :     for (i=l; i<=m; i+=2)
    3420         506 :       a = Fp_mulu(a, i, p);
    3421         344 :     v = Fp_mul(v, k == 1? a: Fp_powu(a, k, p), p);
    3422         344 :     v = gerepileuptoint(av, v);
    3423             :   }
    3424         158 :   return v;
    3425             : }
    3426             : 
    3427             : /*******************************************************************/
    3428             : /**                      LUCAS & FIBONACCI                        **/
    3429             : /*******************************************************************/
    3430             : static void
    3431          56 : lucas(ulong n, GEN *a, GEN *b)
    3432             : {
    3433             :   GEN z, t, zt;
    3434          56 :   if (!n) { *a = gen_2; *b = gen_1; return; }
    3435          49 :   lucas(n >> 1, &z, &t); zt = mulii(z, t);
    3436          49 :   switch(n & 3) {
    3437          14 :     case  0: *a = subiu(sqri(z),2); *b = subiu(zt,1); break;
    3438          14 :     case  1: *a = subiu(zt,1);      *b = addiu(sqri(t),2); break;
    3439           7 :     case  2: *a = addiu(sqri(z),2); *b = addiu(zt,1); break;
    3440          14 :     case  3: *a = addiu(zt,1);      *b = subiu(sqri(t),2);
    3441             :   }
    3442          49 : }
    3443             : 
    3444             : GEN
    3445           7 : fibo(long n)
    3446             : {
    3447           7 :   pari_sp av = avma;
    3448             :   GEN a, b;
    3449           7 :   if (!n) return gen_0;
    3450           7 :   lucas((ulong)(labs(n)-1), &a, &b);
    3451           7 :   a = diviuexact(addii(shifti(a,1),b), 5);
    3452           7 :   if (n < 0 && !odd(n)) setsigne(a, -1);
    3453           7 :   return gerepileuptoint(av, a);
    3454             : }
    3455             : 
    3456             : /*******************************************************************/
    3457             : /*                      CONTINUED FRACTIONS                        */
    3458             : /*******************************************************************/
    3459             : static GEN
    3460     2830744 : icopy_lg(GEN x, long l)
    3461             : {
    3462     2830744 :   long lx = lgefint(x);
    3463             :   GEN y;
    3464             : 
    3465     2830744 :   if (lx >= l) return icopy(x);
    3466          49 :   y = cgeti(l); affii(x, y); return y;
    3467             : }
    3468             : 
    3469             : /* continued fraction of a/b. If y != NULL, stop when partial quotients
    3470             :  * differ from y */
    3471             : static GEN
    3472     2831094 : Qsfcont(GEN a, GEN b, GEN y, ulong k)
    3473             : {
    3474             :   GEN  z, c;
    3475     2831094 :   ulong i, l, ly = lgefint(b);
    3476             : 
    3477             :   /* times 1 / log2( (1+sqrt(5)) / 2 )  */
    3478     2831094 :   l = (ulong)(3 + bit_accuracy_mul(ly, 1.44042009041256));
    3479     2831094 :   if (k > 0 && k+1 > 0 && l > k+1) l = k+1; /* beware overflow */
    3480     2831094 :   if (l > LGBITS) l = LGBITS;
    3481             : 
    3482     2831094 :   z = cgetg(l,t_VEC);
    3483     2831094 :   l--;
    3484     2831094 :   if (y) {
    3485         350 :     pari_sp av = avma;
    3486         350 :     if (l >= (ulong)lg(y)) l = lg(y)-1;
    3487       25209 :     for (i = 1; i <= l; i++)
    3488             :     {
    3489       24985 :       GEN q = gel(y,i);
    3490       24985 :       gel(z,i) = q;
    3491       24985 :       c = b; if (!gequal1(q)) c = mulii(q, b);
    3492       24985 :       c = subii(a, c);
    3493       24985 :       if (signe(c) < 0)
    3494             :       { /* partial quotient too large */
    3495          96 :         c = addii(c, b);
    3496          96 :         if (signe(c) >= 0) i++; /* by 1 */
    3497          96 :         break;
    3498             :       }
    3499       24889 :       if (cmpii(c, b) >= 0)
    3500             :       { /* partial quotient too small */
    3501          30 :         c = subii(c, b);
    3502          30 :         if (cmpii(c, b) < 0) {
    3503             :           /* by 1. If next quotient is 1 in y, add 1 */
    3504          12 :           if (i < l && equali1(gel(y,i+1))) gel(z,i) = addiu(q,1);
    3505          12 :           i++;
    3506             :         }
    3507          30 :         break;
    3508             :       }
    3509       24859 :       if ((i & 0xff) == 0) gerepileall(av, 2, &b, &c);
    3510       24859 :       a = b; b = c;
    3511             :     }
    3512             :   } else {
    3513     2830744 :     a = icopy_lg(a, ly);
    3514     2830744 :     b = icopy(b);
    3515    23447857 :     for (i = 1; i <= l; i++)
    3516             :     {
    3517    23447539 :       gel(z,i) = truedvmdii(a,b,&c);
    3518    23447539 :       if (c == gen_0) { i++; break; }
    3519    20617113 :       affii(c, a); cgiv(c); c = a;
    3520    20617113 :       a = b; b = c;
    3521             :     }
    3522             :   }
    3523     2831094 :   i--;
    3524     2831094 :   if (i > 1 && gequal1(gel(z,i)))
    3525             :   {
    3526         101 :     cgiv(gel(z,i)); --i;
    3527         101 :     gel(z,i) = addui(1, gel(z,i)); /* unclean: leave old z[i] on stack */
    3528             :   }
    3529     2831094 :   setlg(z,i+1); return z;
    3530             : }
    3531             : 
    3532             : static GEN
    3533           0 : sersfcont(GEN a, GEN b, long k)
    3534             : {
    3535           0 :   long i, l = typ(a) == t_POL? lg(a): 3;
    3536             :   GEN y, c;
    3537           0 :   if (lg(b) > l) l = lg(b);
    3538           0 :   if (k > 0 && l > k+1) l = k+1;
    3539           0 :   y = cgetg(l,t_VEC);
    3540           0 :   for (i=1; i<l; i++)
    3541             :   {
    3542           0 :     gel(y,i) = poldivrem(a,b,&c);
    3543           0 :     if (gequal0(c)) { i++; break; }
    3544           0 :     a = b; b = c;
    3545             :   }
    3546           0 :   setlg(y, i); return y;
    3547             : }
    3548             : 
    3549             : GEN
    3550     2831759 : gboundcf(GEN x, long k)
    3551             : {
    3552             :   pari_sp av;
    3553     2831759 :   long tx = typ(x), e;
    3554             :   GEN y, a, b, c;
    3555             : 
    3556     2831759 :   if (k < 0) pari_err_DOMAIN("gboundcf","nmax","<",gen_0,stoi(k));
    3557     2831752 :   if (is_scalar_t(tx))
    3558             :   {
    3559     2831752 :     if (gequal0(x)) return mkvec(gen_0);
    3560     2831647 :     switch(tx)
    3561             :     {
    3562         896 :       case t_INT: return mkveccopy(x);
    3563         357 :       case t_REAL:
    3564         357 :         av = avma;
    3565         357 :         c = mantissa_real(x,&e);
    3566         357 :         if (e < 0) pari_err_PREC("gboundcf");
    3567         350 :         y = int2n(e);
    3568         350 :         a = Qsfcont(c,y, NULL, k);
    3569         350 :         b = addsi(signe(x), c);
    3570         350 :         return gerepilecopy(av, Qsfcont(b,y, a, k));
    3571             : 
    3572     2830394 :       case t_FRAC:
    3573     2830394 :         av = avma;
    3574     2830394 :         return gerepileupto(av, Qsfcont(gel(x,1),gel(x,2), NULL, k));
    3575             :     }
    3576           0 :     pari_err_TYPE("gboundcf",x);
    3577             :   }
    3578             : 
    3579           0 :   switch(tx)
    3580             :   {
    3581           0 :     case t_POL: return mkveccopy(x);
    3582           0 :     case t_SER:
    3583           0 :       av = avma;
    3584           0 :       return gerepileupto(av, gboundcf(ser2rfrac_i(x), k));
    3585           0 :     case t_RFRAC:
    3586           0 :       av = avma;
    3587           0 :       return gerepilecopy(av, sersfcont(gel(x,1), gel(x,2), k));
    3588             :   }
    3589           0 :   pari_err_TYPE("gboundcf",x);
    3590             :   return NULL; /* LCOV_EXCL_LINE */
    3591             : }
    3592             : 
    3593             : static GEN
    3594          14 : sfcont2(GEN b, GEN x, long k)
    3595             : {
    3596          14 :   pari_sp av = avma;
    3597          14 :   long lb = lg(b), tx = typ(x), i;
    3598             :   GEN y,p1;
    3599             : 
    3600          14 :   if (k)
    3601             :   {
    3602           7 :     if (k >= lb) pari_err_DIM("contfrac [too few denominators]");
    3603           0 :     lb = k+1;
    3604             :   }
    3605           7 :   y = cgetg(lb,t_VEC);
    3606           7 :   if (lb==1) return y;
    3607           7 :   if (is_scalar_t(tx))
    3608             :   {
    3609           7 :     if (!is_intreal_t(tx) && tx != t_FRAC) pari_err_TYPE("sfcont2",x);
    3610             :   }
    3611           0 :   else if (tx == t_SER) x = ser2rfrac_i(x);
    3612             : 
    3613           7 :   if (!gequal1(gel(b,1))) x = gmul(gel(b,1),x);
    3614           7 :   for (i = 1;;)
    3615             :   {
    3616          35 :     if (tx == t_REAL)
    3617             :     {
    3618          35 :       long e = expo(x);
    3619          35 :       if (e > 0 && nbits2prec(e+1) > realprec(x)) break;
    3620          35 :       gel(y,i) = floorr(x);
    3621          35 :       p1 = subri(x, gel(y,i));
    3622             :     }
    3623             :     else
    3624             :     {
    3625           0 :       gel(y,i) = gfloor(x);
    3626           0 :       p1 = gsub(x, gel(y,i));
    3627             :     }
    3628          35 :     if (++i >= lb) break;
    3629          28 :     if (gequal0(p1)) break;
    3630          28 :     x = gdiv(gel(b,i),p1);
    3631             :   }
    3632           7 :   setlg(y,i);
    3633           7 :   return gerepilecopy(av,y);
    3634             : }
    3635             : 
    3636             : GEN
    3637         112 : gcf(GEN x) { return gboundcf(x,0); }
    3638             : GEN
    3639           0 : gcf2(GEN b, GEN x) { return contfrac0(x,b,0); }
    3640             : GEN
    3641          49 : contfrac0(GEN x, GEN b, long nmax)
    3642             : {
    3643             :   long tb;
    3644             : 
    3645          49 :   if (!b) return gboundcf(x,nmax);
    3646          28 :   tb = typ(b);
    3647          28 :   if (tb == t_INT) return gboundcf(x,itos(b));
    3648          21 :   if (! is_vec_t(tb)) pari_err_TYPE("contfrac0",b);
    3649          21 :   if (nmax < 0) pari_err_DOMAIN("contfrac","nmax","<",gen_0,stoi(nmax));
    3650          14 :   return sfcont2(b,x,nmax);
    3651             : }
    3652             : 
    3653             : GEN
    3654         252 : contfracpnqn(GEN x, long n)
    3655             : {
    3656         252 :   pari_sp av = avma;
    3657         252 :   long i, lx = lg(x);
    3658             :   GEN M,A,B, p0,p1, q0,q1;
    3659             : 
    3660         252 :   if (lx == 1)
    3661             :   {
    3662          28 :     if (! is_matvec_t(typ(x))) pari_err_TYPE("pnqn",x);
    3663          21 :     if (n >= 0) return cgetg(1,t_MAT);
    3664           7 :     return matid(2);
    3665             :   }
    3666         224 :   switch(typ(x))
    3667             :   {
    3668         182 :     case t_VEC: case t_COL: A = x; B = NULL; break;
    3669          42 :     case t_MAT:
    3670          42 :       switch(lgcols(x))
    3671             :       {
    3672           0 :         case 2: A = row(x,1); B = NULL; break;
    3673          35 :         case 3: A = row(x,2); B = row(x,1); break;
    3674           7 :         default: pari_err_DIM("pnqn [ nbrows != 1,2 ]");
    3675             :                  return NULL; /*LCOV_EXCL_LINE*/
    3676             :       }
    3677          35 :       break;
    3678           0 :     default: pari_err_TYPE("pnqn",x);
    3679             :       return NULL; /*LCOV_EXCL_LINE*/
    3680             :   }
    3681         217 :   p1 = gel(A,1);
    3682         217 :   q1 = B? gel(B,1): gen_1; /* p[0], q[0] */
    3683         217 :   if (n >= 0)
    3684             :   {
    3685         182 :     lx = minss(lx, n+2);
    3686         182 :     if (lx == 2) return gerepilecopy(av, mkmat(mkcol2(p1,q1)));
    3687             :   }
    3688          35 :   else if (lx == 2)
    3689           7 :     return gerepilecopy(av, mkmat2(mkcol2(p1,q1), mkcol2(gen_1,gen_0)));
    3690             :   /* lx >= 3 */
    3691         119 :   p0 = gen_1;
    3692         119 :   q0 = gen_0; /* p[-1], q[-1] */
    3693         119 :   M = cgetg(lx, t_MAT);
    3694         119 :   gel(M,1) = mkcol2(p1,q1);
    3695         399 :   for (i=2; i<lx; i++)
    3696             :   {
    3697         280 :     GEN a = gel(A,i), p2,q2;
    3698         280 :     if (B) {
    3699          84 :       GEN b = gel(B,i);
    3700          84 :       p0 = gmul(b,p0);
    3701          84 :       q0 = gmul(b,q0);
    3702             :     }
    3703         280 :     p2 = gadd(gmul(a,p1),p0); p0=p1; p1=p2;
    3704         280 :     q2 = gadd(gmul(a,q1),q0); q0=q1; q1=q2;
    3705         280 :     gel(M,i) = mkcol2(p1,q1);
    3706             :   }
    3707         119 :   if (n < 0) M = mkmat2(gel(M,lx-1), gel(M,lx-2));
    3708         119 :   return gerepilecopy(av, M);
    3709             : }
    3710             : GEN
    3711           0 : pnqn(GEN x) { return contfracpnqn(x,-1); }
    3712             : /* x = [a0, ..., an] from gboundcf, n >= 0;
    3713             :  * return [[p0, ..., pn], [q0,...,qn]] */
    3714             : GEN
    3715      609308 : ZV_allpnqn(GEN x)
    3716             : {
    3717      609308 :   long i, lx = lg(x);
    3718      609308 :   GEN p0, p1, q0, q1, p2, q2, P,Q, v = cgetg(3,t_VEC);
    3719             : 
    3720      609308 :   gel(v,1) = P = cgetg(lx, t_VEC);
    3721      609308 :   gel(v,2) = Q = cgetg(lx, t_VEC);
    3722      609308 :   p0 = gen_1; q0 = gen_0;
    3723      609308 :   gel(P, 1) = p1 = gel(x,1); gel(Q, 1) = q1 = gen_1;
    3724     2092209 :   for (i=2; i<lx; i++)
    3725             :   {
    3726     1482901 :     GEN a = gel(x,i);
    3727     1482901 :     gel(P, i) = p2 = addmulii(p0, a, p1); p0 = p1; p1 = p2;
    3728     1482901 :     gel(Q, i) = q2 = addmulii(q0, a, q1); q0 = q1; q1 = q2;
    3729             :   }
    3730      609308 :   return v;
    3731             : }
    3732             : 
    3733             : /* write Mod(x,N) as a/b, gcd(a,b) = 1, b <= B (no condition if B = NULL) */
    3734             : static GEN
    3735          42 : mod_to_frac(GEN x, GEN N, GEN B)
    3736             : {
    3737             :   GEN a, b, A;
    3738          42 :   if (B) A = divii(shifti(N, -1), B);
    3739             :   else
    3740             :   {
    3741          14 :     A = sqrti(shifti(N, -1));
    3742          14 :     B = A;
    3743             :   }
    3744          42 :   if (!Fp_ratlift(x, N, A,B,&a,&b) || !equali1( gcdii(a,b) )) return NULL;
    3745          28 :   return equali1(b)? a: mkfrac(a,b);
    3746             : }
    3747             : 
    3748             : static GEN
    3749          70 : mod_to_rfrac(GEN x, GEN N, long B)
    3750             : {
    3751             :   GEN a, b;
    3752          70 :   long A, d = degpol(N);
    3753          70 :   if (B >= 0) A = d-1 - B;
    3754             :   else
    3755             :   {
    3756          42 :     B = d >> 1;
    3757          42 :     A = odd(d)? B : B-1;
    3758             :   }
    3759          70 :   if (varn(N) != varn(x)) x = scalarpol(x, varn(N));
    3760          70 :   if (!RgXQ_ratlift(x, N, A, B, &a,&b) || degpol(RgX_gcd(a,b)) > 0) return NULL;
    3761          56 :   return gdiv(a,b);
    3762             : }
    3763             : 
    3764             : /* k > 0 t_INT, x a t_FRAC, returns the convergent a/b
    3765             :  * of the continued fraction of x with b <= k maximal */
    3766             : static GEN
    3767           7 : bestappr_frac(GEN x, GEN k)
    3768             : {
    3769             :   pari_sp av;
    3770             :   GEN p0, p1, p, q0, q1, q, a, y;
    3771             : 
    3772           7 :   if (cmpii(gel(x,2),k) <= 0) return gcopy(x);
    3773           0 :   av = avma; y = x;
    3774           0 :   p1 = gen_1; p0 = truedvmdii(gel(x,1), gel(x,2), &a); /* = floor(x) */
    3775           0 :   q1 = gen_0; q0 = gen_1;
    3776           0 :   x = mkfrac(a, gel(x,2)); /* = frac(x); now 0<= x < 1 */
    3777             :   for(;;)
    3778             :   {
    3779           0 :     x = ginv(x); /* > 1 */
    3780           0 :     a = typ(x)==t_INT? x: divii(gel(x,1), gel(x,2));
    3781           0 :     if (cmpii(a,k) > 0)
    3782             :     { /* next partial quotient will overflow limits */
    3783             :       GEN n, d;
    3784           0 :       a = divii(subii(k, q1), q0);
    3785           0 :       p = addii(mulii(a,p0), p1); p1=p0; p0=p;
    3786           0 :       q = addii(mulii(a,q0), q1); q1=q0; q0=q;
    3787             :       /* compare |y-p0/q0|, |y-p1/q1| */
    3788           0 :       n = gel(y,1);
    3789           0 :       d = gel(y,2);
    3790           0 :       if (abscmpii(mulii(q1, subii(mulii(q0,n), mulii(d,p0))),
    3791             :                    mulii(q0, subii(mulii(q1,n), mulii(d,p1)))) < 0)
    3792           0 :                    { p1 = p0; q1 = q0; }
    3793           0 :       break;
    3794             :     }
    3795           0 :     p = addii(mulii(a,p0), p1); p1=p0; p0=p;
    3796           0 :     q = addii(mulii(a,q0), q1); q1=q0; q0=q;
    3797             : 
    3798           0 :     if (cmpii(q0,k) > 0) break;
    3799           0 :     x = gsub(x,a); /* 0 <= x < 1 */
    3800           0 :     if (typ(x) == t_INT) { p1 = p0; q1 = q0; break; } /* x = 0 */
    3801             : 
    3802             :   }
    3803           0 :   return gerepileupto(av, gdiv(p1,q1));
    3804             : }
    3805             : /* k > 0 t_INT, x != 0 a t_REAL, returns the convergent a/b
    3806             :  * of the continued fraction of x with b <= k maximal */
    3807             : static GEN
    3808     1227891 : bestappr_real(GEN x, GEN k)
    3809             : {
    3810     1227891 :   pari_sp av = avma;
    3811     1227891 :   GEN kr, p0, p1, p, q0, q1, q, a, y = x;
    3812             : 
    3813     1227891 :   p1 = gen_1; a = p0 = floorr(x);
    3814     1227812 :   q1 = gen_0; q0 = gen_1;
    3815     1227812 :   x = subri(x,a); /* 0 <= x < 1 */
    3816     1227854 :   if (!signe(x)) { cgiv(x); return a; }
    3817     1114133 :   kr = itor(k, realprec(x));
    3818             :   for(;;)
    3819     1196075 :   {
    3820             :     long d;
    3821     2310230 :     x = invr(x); /* > 1 */
    3822     2309960 :     if (cmprr(x,kr) > 0)
    3823             :     { /* next partial quotient will overflow limits */
    3824     1099110 :       a = divii(subii(k, q1), q0);
    3825     1099107 :       p = addii(mulii(a,p0), p1); p1=p0; p0=p;
    3826     1099154 :       q = addii(mulii(a,q0), q1); q1=q0; q0=q;
    3827             :       /* compare |y-p0/q0|, |y-p1/q1| */
    3828     1099130 :       if (abscmprr(mulir(q1, subri(mulir(q0,y), p0)),
    3829             :                    mulir(q0, subri(mulir(q1,y), p1))) < 0)
    3830      131215 :                    { p1 = p0; q1 = q0; }
    3831     1099188 :       break;
    3832             :     }
    3833     1210981 :     d = nbits2prec(expo(x) + 1);
    3834     1210981 :     if (d > realprec(x)) { p1 = p0; q1 = q0; break; } /* original x was ~ 0 */
    3835             : 
    3836     1210791 :     a = truncr(x); /* truncr(x) will NOT raise e_PREC */
    3837     1210760 :     p = addii(mulii(a,p0), p1); p1=p0; p0=p;
    3838     1210772 :     q = addii(mulii(a,q0), q1); q1=q0; q0=q;
    3839             : 
    3840     1210780 :     if (cmpii(q0,k) > 0) break;
    3841     1204560 :     x = subri(x,a); /* 0 <= x < 1 */
    3842     1204572 :     if (!signe(x)) { p1 = p0; q1 = q0; break; }
    3843             :   }
    3844     1114092 :   if (signe(q1) < 0) { togglesign_safe(&p1); togglesign_safe(&q1); }
    3845     1114092 :   return gerepilecopy(av, equali1(q1)? p1: mkfrac(p1,q1));
    3846             : }
    3847             : 
    3848             : /* k t_INT or NULL */
    3849             : static GEN
    3850     2142646 : bestappr_Q(GEN x, GEN k)
    3851             : {
    3852     2142646 :   long lx, tx = typ(x), i;
    3853             :   GEN a, y;
    3854             : 
    3855     2142646 :   switch(tx)
    3856             :   {
    3857         154 :     case t_INT: return icopy(x);
    3858           7 :     case t_FRAC: return k? bestappr_frac(x, k): gcopy(x);
    3859     1460464 :     case t_REAL:
    3860     1460464 :       if (!signe(x)) return gen_0;
    3861             :       /* i <= e iff nbits2lg(e+1) > lg(x) iff floorr(x) fails */
    3862     1227878 :       i = bit_prec(x); if (i <= expo(x)) return NULL;
    3863     1227891 :       return bestappr_real(x, k? k: int2n(i));
    3864             : 
    3865          28 :     case t_INTMOD: {
    3866          28 :       pari_sp av = avma;
    3867          28 :       a = mod_to_frac(gel(x,2), gel(x,1), k); if (!a) return NULL;
    3868          21 :       return gerepilecopy(av, a);
    3869             :     }
    3870          14 :     case t_PADIC: {
    3871          14 :       pari_sp av = avma;
    3872          14 :       long v = valp(x);
    3873          14 :       a = mod_to_frac(gel(x,4), gel(x,3), k); if (!a) return NULL;
    3874           7 :       if (v) a = gmul(a, powis(gel(x,2), v));
    3875           7 :       return gerepilecopy(av, a);
    3876             :     }
    3877             : 
    3878         294 :     case t_COMPLEX: {
    3879         294 :       pari_sp av = avma;
    3880         294 :       y = cgetg(3, t_COMPLEX);
    3881         294 :       gel(y,2) = bestappr(gel(x,2), k);
    3882         294 :       gel(y,1) = bestappr(gel(x,1), k);
    3883         294 :       if (gequal0(gel(y,2))) return gerepileupto(av, gel(y,1));
    3884          91 :       return y;
    3885             :     }
    3886           0 :     case t_SER:
    3887           0 :       if (ser_isexactzero(x)) return gcopy(x);
    3888             :       /* fall through */
    3889             :     case t_POLMOD: case t_POL: case t_RFRAC:
    3890             :     case t_VEC: case t_COL: case t_MAT:
    3891      681685 :       y = cgetg_copy(x, &lx);
    3892      681725 :       if (lontyp[tx] == 1) i = 1; else { y[1] = x[1]; i = 2; }
    3893     2757529 :       for (; i<lx; i++)
    3894             :       {
    3895     2075829 :         a = bestappr_Q(gel(x,i),k); if (!a) return NULL;
    3896     2075804 :         gel(y,i) = a;
    3897             :       }
    3898      681700 :       if (tx == t_POL) return normalizepol(y);
    3899      681686 :       if (tx == t_SER) return normalizeser(y);
    3900      681686 :       return y;
    3901             :   }
    3902           0 :   pari_err_TYPE("bestappr_Q",x);
    3903             :   return NULL; /* LCOV_EXCL_LINE */
    3904             : }
    3905             : 
    3906             : static GEN
    3907          56 : bestappr_ser(GEN x, long B)
    3908             : {
    3909          56 :   long dN, v = valser(x), lx = lg(x);
    3910             :   GEN t;
    3911          56 :   x = normalizepol(ser2pol_i(x, lx));
    3912          56 :   dN = lx-2;
    3913          56 :   if (v > 0)
    3914             :   {
    3915          14 :     x = RgX_shift_shallow(x, v);
    3916          14 :     dN += v;
    3917             :   }
    3918          42 :   else if (v < 0)
    3919             :   {
    3920           7 :     if (B >= 0) B = maxss(B+v, 0);
    3921             :   }
    3922          56 :   t = mod_to_rfrac(x, pol_xn(dN, varn(x)), B);
    3923          56 :   if (!t) return NULL;
    3924          42 :   if (v < 0)
    3925             :   {
    3926             :     GEN a, b;
    3927             :     long vx;
    3928           7 :     if (typ(t) == t_POL) return RgX_mulXn(t, v);
    3929             :     /* t_RFRAC */
    3930           7 :     vx = varn(x);
    3931           7 :     a = gel(t,1);
    3932           7 :     b = gel(t,2);
    3933           7 :     v -= RgX_valrem(b, &b);
    3934           7 :     if (typ(a) == t_POL && varn(a) == vx) v += RgX_valrem(a, &a);
    3935           7 :     if (v < 0) b = RgX_shift(b, -v);
    3936           0 :     else if (v > 0) {
    3937           0 :       if (typ(a) != t_POL || varn(a) != vx) a = scalarpol_shallow(a, vx);
    3938           0 :       a = RgX_shift(a, v);
    3939             :     }
    3940           7 :     t = mkrfraccopy(a, b);
    3941             :   }
    3942          42 :   return t;
    3943             : }
    3944             : static GEN bestappr_RgX(GEN x, long B);
    3945             : /* x t_POLMOD, B >= 0 or < 0 [omit condition on B].
    3946             :  * Look for coprime t_POL a,b, deg(b)<=B, such that a/b = x */
    3947             : static GEN
    3948          77 : bestappr_RgX(GEN x, long B)
    3949             : {
    3950          77 :   long i, lx, tx = typ(x);
    3951             :   GEN y, t;
    3952          77 :   switch(tx)
    3953             :   {
    3954           0 :     case t_INT: case t_REAL: case t_INTMOD: case t_FRAC:
    3955             :     case t_COMPLEX: case t_PADIC: case t_QUAD: case t_POL:
    3956           0 :       return gcopy(x);
    3957             : 
    3958          14 :     case t_RFRAC: {
    3959          14 :       pari_sp av = avma;
    3960          14 :       if (B < 0 || degpol(gel(x,2)) <= B) return gcopy(x);
    3961           7 :       x = rfrac_to_ser_i(x, 2*B+1);
    3962           7 :       t = bestappr_ser(x, B); if (!t) return NULL;
    3963           0 :       return gerepileupto(av, t);
    3964             :     }
    3965          14 :     case t_POLMOD: {
    3966          14 :       pari_sp av = avma;
    3967          14 :       t = mod_to_rfrac(gel(x,2), gel(x,1), B); if (!t) return NULL;
    3968          14 :       return gerepileupto(av, t);
    3969             :     }
    3970          49 :     case t_SER: {
    3971          49 :       pari_sp av = avma;
    3972          49 :       t = bestappr_ser(x, B); if (!t) return NULL;
    3973          42 :       return gerepileupto(av, t);
    3974             :     }
    3975             : 
    3976           0 :     case t_VEC: case t_COL: case t_MAT:
    3977           0 :       y = cgetg_copy(x, &lx);
    3978           0 :       if (lontyp[tx] == 1) i = 1; else { y[1] = x[1]; i = 2; }
    3979           0 :       for (; i<lx; i++)
    3980             :       {
    3981           0 :         t = bestappr_RgX(gel(x,i),B); if (!t) return NULL;
    3982           0 :         gel(y,i) = t;
    3983             :       }
    3984           0 :       return y;
    3985             :   }
    3986           0 :   pari_err_TYPE("bestappr_RgX",x);
    3987             :   return NULL; /* LCOV_EXCL_LINE */
    3988             : }
    3989             : 
    3990             : /* allow k = NULL: maximal accuracy */
    3991             : GEN
    3992       66805 : bestappr(GEN x, GEN k)
    3993             : {
    3994       66805 :   pari_sp av = avma;
    3995       66805 :   if (k) { /* replace by floor(k) */
    3996       66483 :     switch(typ(k))
    3997             :     {
    3998        4368 :       case t_INT:
    3999        4368 :         break;
    4000       62115 :       case t_REAL: case t_FRAC:
    4001       62115 :         k = floor_safe(k); /* left on stack for efficiency */
    4002       62116 :         if (!signe(k)) k = gen_1;
    4003       62116 :         break;
    4004           0 :       default:
    4005           0 :         pari_err_TYPE("bestappr [bound type]", k);
    4006           0 :         break;
    4007             :     }
    4008         322 :   }
    4009       66806 :   x = bestappr_Q(x, k);
    4010       66806 :   if (!x) { set_avma(av); return cgetg(1,t_VEC); }
    4011       66792 :   return x;
    4012             : }
    4013             : GEN
    4014          77 : bestapprPade(GEN x, long B)
    4015             : {
    4016          77 :   pari_sp av = avma;
    4017          77 :   GEN t = bestappr_RgX(x, B);
    4018          77 :   if (!t) { set_avma(av); return cgetg(1,t_VEC); }
    4019          63 :   return t;
    4020             : }
    4021             : 
    4022             : GEN
    4023          42 : serPade(GEN S, long p, long q)
    4024             : {
    4025          42 :   pari_sp ltop = avma;
    4026          42 :   GEN M, B, Q, P, y, LT = gen_1;
    4027          42 :   long c, va, tS = typ(S), i, j, lS, val;
    4028          42 :   if (p < 0 || q < 0) pari_err(e_MISC, "need p and q nonnegative in serPade");
    4029          42 :   switch(tS)
    4030             :   {
    4031           0 :     case t_INT: case t_REAL: case t_INTMOD: case t_FRAC:
    4032             :     case t_COMPLEX: case t_PADIC: case t_QUAD: case t_POLMOD:
    4033           0 :       return gcopy(S);
    4034           7 :     case t_RFRAC: S = rfrac_to_ser_i(S, p + q + 4); break;
    4035          14 :     case t_POL: S = RgX_to_ser(S, p + q + 4); break;
    4036          21 :     case t_SER: break;
    4037           0 :     case t_VEC: case t_COL: case t_MAT:
    4038           0 :       y = cgetg_copy(S, &lS);
    4039           0 :       if (lontyp[tS] == 1) i = 1; else { y[1] = S[1]; i = 2; }
    4040           0 :       for (; i<lS; i++) gel(y,i) = serPade(gel(S,i),p,q);
    4041           0 :       return y;
    4042           0 :     default: pari_err_TYPE("serPade", S);
    4043             :   }
    4044          42 :   va = gvar(S);
    4045          42 :   val = valser(S);
    4046          42 :   if (val) { LT = gpowgs(pol_x(va), val); S = gdiv(S, LT); }
    4047          42 :   if (val < 0)
    4048             :   {
    4049           7 :     if (q < -val) pari_err(e_MISC, "-valuation larger than q in serPade");
    4050           7 :     q += val;
    4051             :   }
    4052          42 :   if (val > 0)
    4053             :   {
    4054           7 :     if (p < val) pari_err(e_MISC, "+valuation larger than p in serPade");
    4055           7 :     p -= val;
    4056             :   }
    4057          42 :   if (p + q >= serprec(S, va))
    4058           0 :     pari_err(e_MISC, "insufficient series precision in serPade");
    4059          42 :   M = cgetg(q + 2, t_MAT);
    4060         217 :   for (j = 1; j <= q + 1; j++)
    4061             :   {
    4062         175 :     GEN C = cgetg(q + 1, t_COL);
    4063         749 :     for (i = 1; i <= q; i++) gel(C, i) = polcoeff0(S, i - j + p + 1, va);
    4064         175 :     gel(M, j) = C;
    4065             :   }
    4066          42 :   B = gel(ker(M), 1);
    4067          49 :   c = 1; while (gequal0(gel(B, c))) c++; B = gdiv(B, gel(B, c));
    4068          42 :   Q = gtopolyrev(B, va);
    4069          42 :   M = cgetg(q + 2, t_MAT);
    4070         217 :   for (j = 1; j <= q + 1; j++)
    4071             :   {
    4072         175 :     GEN C = cgetg(p + 2, t_COL);
    4073         945 :     for (i = 1; i <= p + 1; i++) gel(C, i) = polcoeff0(S, i - j, va);
    4074         175 :     gel(M, j) = C;
    4075             :   }
    4076          42 :   P = gtopolyrev(gmul(M, B), va);
    4077          42 :   return gerepileupto(ltop, val ? gmul(LT, gdiv(P, Q)) : gdiv(P, Q));
    4078             : }

Generated by: LCOV version 1.14