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 29419-8afb0ed749) Lines: 2072 2272 91.2 %
Date: 2024-07-02 09:03:41 Functions: 216 230 93.9 %
Legend: Lines: hit not hit

          Line data    Source code
       1             : /* Copyright (C) 2000  The PARI group.
       2             : 
       3             : This file is part of the PARI/GP package.
       4             : 
       5             : PARI/GP is free software; you can redistribute it and/or modify it under the
       6             : terms of the GNU General Public License as published by the Free Software
       7             : Foundation; either version 2 of the License, or (at your option) any later
       8             : version. It is distributed in the hope that it will be useful, but WITHOUT
       9             : ANY WARRANTY WHATSOEVER.
      10             : 
      11             : Check the License for details. You should have received a copy of it, along
      12             : with the package; see the file 'COPYING'. If not, write to the Free Software
      13             : Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. */
      14             : 
      15             : /*********************************************************************/
      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      464970 : 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      464970 : 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      532926 : u_is_gener_expo(ulong p, GEN L0)
      54             : {
      55      532926 :   const ulong q = p >> 1;
      56             :   long i;
      57             :   GEN L;
      58      532926 :   if (!L0) L0 = u_odd_prime_divisors(q);
      59      532928 :   L = cgetg_copy(L0,&i);
      60     1155480 :   while (--i) L[i] = q / uel(L0,i);
      61      532927 :   return L;
      62             : }
      63             : 
      64             : int
      65     1629979 : is_gener_Fl(ulong x, ulong p, ulong p_1, GEN L)
      66             : {
      67             :   long i;
      68     1629979 :   if (krouu(x, p) >= 0) return 0;
      69     1380579 :   for (i=lg(L)-1; i; i--)
      70             :   {
      71      842804 :     ulong t = Fl_powu(x, uel(L,i), p);
      72      842804 :     if (t == p_1 || t == 1) return 0;
      73             :   }
      74      537775 :   return 1;
      75             : }
      76             : /* assume p prime */
      77             : ulong
      78     1056222 : pgener_Fl_local(ulong p, GEN L0)
      79             : {
      80     1056222 :   const pari_sp av = avma;
      81     1056222 :   const ulong p_1 = p-1;
      82             :   long x;
      83             :   GEN L;
      84     1056222 :   if (p <= 19) switch(p)
      85             :   { /* quick trivial cases */
      86          28 :     case 2:  return 1;
      87      116618 :     case 7:
      88      116618 :     case 17: return 3;
      89      406694 :     default: return 2;
      90             :   }
      91      532882 :   L = u_is_gener_expo(p,L0);
      92     1622127 :   for (x = 2;; x++)
      93     1622127 :     if (is_gener_Fl(x,p,p_1,L)) return gc_ulong(av, x);
      94             : }
      95             : ulong
      96      581350 : 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       14042 : is_gener_Fp(GEN x, GEN p, GEN p_1, GEN L)
     102             : {
     103       14042 :   long i, t = lgefint(x)==3? kroui(x[2], p): kronecker(x, p);
     104       14042 :   if (t >= 0) return 0;
     105       22102 :   for (i = lg(L)-1; i; i--)
     106             :   {
     107       14571 :     GEN t = Fp_pow(x, gel(L,i), p);
     108       14571 :     if (equalii(t, p_1) || equali1(t)) return 0;
     109             :   }
     110        7531 :   return 1;
     111             : }
     112             : 
     113             : /* assume p prime, return a generator of all L[i]-Sylows in F_p^*. */
     114             : GEN
     115      357911 : pgener_Fp_local(GEN p, GEN L0)
     116             : {
     117      357911 :   pari_sp av0 = avma;
     118             :   GEN x, p_1, L;
     119      357911 :   if (lgefint(p) == 3)
     120             :   {
     121             :     ulong z;
     122      353008 :     if (p[2] == 2) return gen_1;
     123      258069 :     if (L0) L0 = ZV_to_nv(L0);
     124      258067 :     z = pgener_Fl_local(uel(p,2), L0);
     125      258096 :     return gc_utoipos(av0, z);
     126             :   }
     127        4903 :   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       44163 : pgener_Fp(GEN p) { return pgener_Fp_local(p, NULL); }
     135             : 
     136             : ulong
     137      205262 : pgener_Zl(ulong p)
     138             : {
     139      205262 :   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      205262 :   if (p == 40487) return 10;
     142             : #ifndef LONG_IS_64BIT
     143       29744 :   return pgener_Fl(p);
     144             : #else
     145      175518 :   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      170506 : pgener_Zp(GEN p)
     162             : {
     163      170506 :   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      214613 : rootsof1_Fl(ulong n, ulong p)
     251             : {
     252      214613 :   pari_sp av = avma;
     253      214613 :   GEN L = u_odd_prime_divisors(n); /* 2 implicit in pgener_Fl_local */
     254      214613 :   ulong z = pgener_Fl_local(p, L);
     255      214613 :   z = Fl_powu(z, (p-1) / n, p); /* prim. n-th root of 1 */
     256      214613 :   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   320256285 : ome(long t)
     364             : {
     365   320256285 :   switch(t & 7)
     366             :   {
     367   181567689 :     case 3:
     368   181567689 :     case 5: return 1;
     369   138688596 :     default: return 0;
     370             :   }
     371             : }
     372             : /* t a t_INT, is t = 3,5 mod 8 ? */
     373             : static int
     374     5590128 : gome(GEN t)
     375     5590128 : { return signe(t)? ome( mod2BIL(t) ): 0; }
     376             : 
     377             : /* assume y odd, return kronecker(x,y) * s */
     378             : static long
     379   226866817 : krouu_s(ulong x, ulong y, long s)
     380             : {
     381   226866817 :   ulong x1 = x, y1 = y, z;
     382  1028727875 :   while (x1)
     383             :   {
     384   801928464 :     long r = vals(x1);
     385   801914897 :     if (r)
     386             :     {
     387   426408868 :       if (odd(r) && ome(y1)) s = -s;
     388   426355029 :       x1 >>= r;
     389             :     }
     390   801861058 :     if (x1 & y1 & 2) s = -s;
     391   801861058 :     z = y1 % x1; y1 = x1; x1 = z;
     392             :   }
     393   226799411 :   return (y1 == 1)? s: 0;
     394             : }
     395             : 
     396             : long
     397    11943670 : kronecker(GEN x, GEN y)
     398             : {
     399    11943670 :   pari_sp av = avma;
     400    11943670 :   long s = 1, r;
     401             :   ulong xu;
     402             : 
     403    11943670 :   if (typ(x) != t_INT) pari_err_TYPE("kronecker",x);
     404    11943670 :   if (typ(y) != t_INT) pari_err_TYPE("kronecker",y);
     405    11943670 :   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    11943537 :   r = vali(y);
     411    11943538 :   if (r)
     412             :   {
     413     1348898 :     if (!mpodd(x)) return gc_long(av,0);
     414      321711 :     if (odd(r) && gome(x)) s = -s;
     415      321711 :     y = shifti(y,-r);
     416             :   }
     417    10916351 :   x = modii(x,y);
     418    13302267 :   while (lgefint(x) > 3) /* x < y */
     419             :   {
     420             :     GEN z;
     421     2385995 :     r = vali(x);
     422     2385808 :     if (r)
     423             :     {
     424     1303306 :       if (odd(r) && gome(y)) s = -s;
     425     1303255 :       x = shifti(x,-r);
     426             :     }
     427             :     /* x=3 mod 4 && y=3 mod 4 ? (both are odd here) */
     428     2385361 :     if (mod2BIL(x) & mod2BIL(y) & 2) s = -s;
     429     2384845 :     z = remii(y,x); y = x; x = z;
     430     2385903 :     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    10916272 :   xu = itou(x);
     437    10916269 :   if (!xu) return is_pm1(y)? s: 0;
     438    10819071 :   r = vals(xu);
     439    10819065 :   if (r)
     440             :   {
     441     5747122 :     if (odd(r) && gome(y)) s = -s;
     442     5747122 :     xu >>= r;
     443             :   }
     444             :   /* x=3 mod 4 && y=3 mod 4 ? (both are odd here) */
     445    10819065 :   if (xu & mod2BIL(y) & 2) s = -s;
     446    10819056 :   return gc_long(av, krouu_s(umodiu(y,xu), xu, s));
     447             : }
     448             : 
     449             : long
     450       39641 : krois(GEN x, long y)
     451             : {
     452             :   ulong yu;
     453       39641 :   long s = 1;
     454             : 
     455       39641 :   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       39613 :     yu = (ulong)y;
     462       39613 :   if (!odd(yu))
     463             :   {
     464             :     long r;
     465       18361 :     if (!mpodd(x)) return 0;
     466       12467 :     r = vals(yu); yu >>= r;
     467       12467 :     if (odd(r) && gome(x)) s = -s;
     468             :   }
     469       33719 :   return krouu_s(umodiu(x, yu), yu, s);
     470             : }
     471             : /* assume y != 0 */
     472             : long
     473    27524958 : kroiu(GEN x, ulong y)
     474             : {
     475             :   long r;
     476    27524958 :   if (odd(y)) return krouu_s(umodiu(x,y), y, 1);
     477      302152 :   if (!mpodd(x)) return 0;
     478      207485 :   r = vals(y); y >>= r;
     479      207493 :   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      179625 : krouodd(ulong x, GEN y, long s)
     485             : {
     486             :   long r;
     487      179625 :   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       15243 :     if (odd(r) && gome(y)) s = -s;
     493       15243 :     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      144537 : krosi(long x, GEN y)
     502             : {
     503      144537 :   const pari_sp av = avma;
     504      144537 :   long s = 1, r;
     505      144537 :   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      144537 :   r = vali(y);
     511      144537 :   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      144537 :   if (x < 0) { x = -x; if (mod4(y) == 3) s = -s; }
     518      144537 :   return gc_long(av, krouodd((ulong)x, y, s));
     519             : }
     520             : 
     521             : long
     522       35088 : kroui(ulong x, GEN y)
     523             : {
     524       35088 :   const pari_sp av = avma;
     525       35088 :   long s = 1, r;
     526       35088 :   switch (signe(y))
     527             :   {
     528           0 :     case -1: y = negi(y); break;
     529           0 :     case 0: return x==1UL;
     530             :   }
     531       35088 :   r = vali(y);
     532       35088 :   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       35088 :   return gc_long(av, krouodd(x, y, s));
     539             : }
     540             : 
     541             : long
     542    97730389 : kross(long x, long y)
     543             : {
     544             :   ulong yu;
     545    97730389 :   long s = 1;
     546             : 
     547    97730389 :   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    97661446 :     yu = (ulong)y;
     554    97730361 :   if (!odd(yu))
     555             :   {
     556             :     long r;
     557    23570708 :     if (!odd(x)) return 0;
     558    16652168 :     r = vals(yu); yu >>= r;
     559    16652168 :     if (odd(r) && ome(x)) s = -s;
     560             :   }
     561    90811821 :   x %= (long)yu; if (x < 0) x += yu;
     562    90811821 :   return krouu_s((ulong)x, yu, s);
     563             : }
     564             : 
     565             : long
     566    97617223 : krouu(ulong x, ulong y)
     567             : {
     568             :   long r;
     569    97617223 :   if (odd(y)) return krouu_s(x, y, 1);
     570       12793 :   if (!odd(x)) return 0;
     571       12962 :   r = vals(y); y >>= r;
     572       12962 :   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      135373 : hilbertii(GEN x, GEN y, GEN p)
     589             : {
     590             :   pari_sp av;
     591             :   long oddvx, oddvy, z;
     592             : 
     593      135373 :   if (!p) return mphilbertoo(x,y);
     594      128065 :   if (is_pm1(p) || signe(p) < 0) pari_err_PRIME("hilbertii",p);
     595      128065 :   if (!signe(x) || !signe(y)) return 0;
     596      128044 :   av = avma;
     597      128044 :   oddvx = odd(Z_pvalrem(x,p,&x));
     598      128044 :   oddvy = odd(Z_pvalrem(y,p,&y));
     599             :   /* x, y are p-units, compute hilbert(x * p^oddvx, y * p^oddvy, p) */
     600      128044 :   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      115941 :     z = (oddvx && oddvy && mod4(p) == 3)? -1: 1;
     609      115941 :     if (oddvx && kronecker(y,p) < 0) z = -z;
     610      115941 :     if (oddvy && kronecker(x,p) < 0) z = -z;
     611             :   }
     612      128044 :   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       60445 : hilbert(GEN x, GEN y, GEN p)
     664             : {
     665       60445 :   pari_sp av = avma;
     666       60445 :   long tx = typ(x), ty = typ(y);
     667             : 
     668       60445 :   if (p && typ(p) != t_INT) pari_err_TYPE("hilbert",p);
     669       60445 :   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       60368 :   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       60354 :   if (tx == t_INTMOD) { x = lift_intmod(x, &p); tx = t_INT; }
     692       60151 :   if (ty == t_INTMOD) { y = lift_intmod(y, &p); ty = t_INT; }
     693             : 
     694       60095 :   if (tx == t_PADIC) { x = lift_padic(x, &p); tx = t_INT; }
     695       60032 :   if (ty == t_PADIC) { y = lift_padic(y, &p); ty = t_INT; }
     696             : 
     697       59955 :   if (tx == t_FRAC) { tx = t_INT; x = p? mulii(gel(x,1),gel(x,2)): gel(x,1); }
     698       59955 :   if (ty == t_FRAC) { ty = t_INT; y = p? mulii(gel(y,1),gel(y,2)): gel(y,1); }
     699             : 
     700       59955 :   if (tx != t_INT || ty != t_INT) pari_err_TYPE2("hilbert",x,y);
     701       59955 :   if (p && !signe(p)) p = NULL;
     702       59955 :   return gc_long(av, hilbertii(x,y,p));
     703             : }
     704             : 
     705             : /*******************************************************************/
     706             : /*                       SQUARE ROOT MODULO p                      */
     707             : /*******************************************************************/
     708             : static void
     709     2250328 : checkp(ulong q, ulong p)
     710     2250328 : { 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    11588014 : nonsquare1_Fl(ulong p)
     714             : {
     715             :   forprime_t S;
     716             :   ulong q;
     717    11588014 :   if ((p & 7UL) != 1) return 2UL;
     718     4348849 :   q = p % 3; if (q == 2) return 3UL;
     719     1396691 :   checkp(q, p);
     720     1403693 :   q = p % 5; if (q == 2 || q == 3) return 5UL;
     721      523876 :   checkp(q, p);
     722      523887 :   q = p % 7; if (q != 4 && q >= 3) return 7UL;
     723      195804 :   checkp(q, p);
     724             :   /* log^2(2^64) < 1968 is enough under GRH (and p^(1/4)log(p) without it)*/
     725      195843 :   u_forprime_init(&S, 11, 1967);
     726      323063 :   while ((q = u_forprime_next(&S)))
     727             :   {
     728      323059 :     if (krouu(q, p) < 0) return q;
     729      127211 :     checkp(q, p);
     730             :   }
     731           0 :   checkp(0, p);
     732             :   return 0; /*LCOV_EXCL_LINE*/
     733             : }
     734             : /* p > 2 a prime */
     735             : ulong
     736        7935 : nonsquare_Fl(ulong p)
     737        7935 : { return ((p & 3UL) == 3)? p-1: nonsquare1_Fl(p); }
     738             : 
     739             : /* allow pi = 0 */
     740             : ulong
     741      177394 : Fl_2gener_pre(ulong p, ulong pi)
     742             : {
     743      177394 :   ulong p1 = p-1;
     744      177394 :   long e = vals(p1);
     745      177380 :   if (e == 1) return p1;
     746       65811 :   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    12213060 : Fl_sqrt_i(ulong a, ulong y, ulong p)
     760             : {
     761             :   long i, e, k;
     762             :   ulong p1, q, v, w;
     763             : 
     764    12213060 :   if (!a) return 0;
     765    10940685 :   p1 = p - 1; e = vals(p1);
     766    10941304 :   if (e == 0) /* p = 2 */
     767             :   {
     768      649978 :     if (p != 2) pari_err_PRIME("Fl_sqrt [modulus]",utoi(p));
     769      650969 :     return ((a & 1) == 0)? 0: 1;
     770             :   }
     771    10291326 :   if (e == 1)
     772             :   {
     773     4853976 :     v = Fl_powu(a, (p+1) >> 2, p);
     774     4854037 :     if (Fl_sqr(v, p) != a) return ~0UL;
     775     4849156 :     p1 = p - v; if (v > p1) v = p1;
     776     4849156 :     return v;
     777             :   }
     778     5437350 :   q = p1 >> e; /* q = (p-1)/2^oo is odd */
     779     5437350 :   p1 = Fl_powu(a, q >> 1, p); /* a ^ [(q-1)/2] */
     780     5437358 :   if (!p1) return 0;
     781     5437358 :   v = Fl_mul(a, p1, p);
     782     5437354 :   w = Fl_mul(v, p1, p);
     783     5437344 :   if (!y) y = Fl_powu(nonsquare1_Fl(p), q, p);
     784     9242379 :   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     3806915 :     p1 = Fl_sqr(w, p);
     788     6323765 :     for (k=1; p1 != 1 && k < e; k++) p1 = Fl_sqr(p1, p);
     789     3806920 :     if (k == e) return ~0UL;
     790             :     /* w ^ (2^k) = 1 --> w = y ^ (u * 2^(e-k)), u odd */
     791     3804843 :     p1 = y;
     792     5080733 :     for (i=1; i < e-k; i++) p1 = Fl_sqr(p1, p);
     793     3804842 :     y = Fl_sqr(p1, p); e = k;
     794     3804878 :     w = Fl_mul(y, w, p);
     795     3804875 :     v = Fl_mul(v, p1, p);
     796             :   }
     797     5435464 :   p1 = p - v; if (v > p1) v = p1;
     798     5435464 :   return v;
     799             : }
     800             : 
     801             : /* Tonelli-Shanks. Assume p is prime and (a,p) != -1. Allow pi = 0 */
     802             : ulong
     803    33474281 : 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    33474281 :   if (!pi) return Fl_sqrt_i(a, y, p);
     809    21261166 :   if (!a) return 0;
     810    21138285 :   p1 = p - 1; e = vals(p1);
     811    21137274 :   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    21152021 :   if (e == 1)
     817             :   {
     818    15015011 :     v = Fl_powu_pre(a, (p+1) >> 2, p, pi);
     819    15001679 :     if (Fl_sqr_pre(v, p, pi) != a) return ~0UL;
     820    15007017 :     p1 = p - v; if (v > p1) v = p1;
     821    15007017 :     return v;
     822             :   }
     823     6137010 :   q = p1 >> e; /* q = (p-1)/2^oo is odd */
     824     6137010 :   p1 = Fl_powu_pre(a, q >> 1, p, pi); /* a ^ [(q-1)/2] */
     825     6133805 :   if (!p1) return 0;
     826     6133805 :   v = Fl_mul_pre(a, p1, p, pi);
     827     6134536 :   w = Fl_mul_pre(v, p1, p, pi);
     828     6134160 :   if (!y) y = Fl_powu_pre(nonsquare1_Fl(p), q, p, pi);
     829    11669476 :   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     5534319 :     p1 = Fl_sqr_pre(w,p,pi);
     833    10345937 :     for (k=1; p1 != 1 && k < e; k++) p1 = Fl_sqr_pre(p1,p,pi);
     834     5533423 :     if (k == e) return ~0UL;
     835             :     /* w ^ (2^k) = 1 --> w = y ^ (u * 2^(e-k)), u odd */
     836     5533331 :     p1 = y;
     837     7267676 :     for (i=1; i < e-k; i++) p1 = Fl_sqr_pre(p1, p, pi);
     838     5533484 :     y = Fl_sqr_pre(p1, p, pi); e = k;
     839     5536177 :     w = Fl_mul_pre(y, w, p, pi);
     840     5534559 :     v = Fl_mul_pre(v, p1, p, pi);
     841             :   }
     842     6135157 :   p1 = p - v; if (v > p1) v = p1;
     843     6135157 :   return v;
     844             : }
     845             : 
     846             : ulong
     847    12269035 : Fl_sqrt(ulong a, ulong p)
     848    12269035 : { ulong pi = (p & HIGHMASK)? get_Fl_red(p): 0; return Fl_sqrt_pre_i(a, 0, p, pi); }
     849             : 
     850             : ulong
     851    21044754 : Fl_sqrt_pre(ulong a, ulong p, ulong pi)
     852    21044754 : { return Fl_sqrt_pre_i(a, 0, p, pi); }
     853             : 
     854             : /* allow pi = 0 */
     855             : static ulong
     856      141245 : Fl_lgener_pre_all(ulong l, long e, ulong r, ulong p, ulong pi, ulong *pt_m)
     857             : {
     858      141245 :   ulong x, y, m, le1 = upowuu(l, e-1);
     859      141245 :   for (x = 2; ; x++)
     860             :   {
     861      172098 :     y = Fl_powu_pre(x, r, p, pi);
     862      172095 :     if (y==1) continue;
     863      154980 :     m = Fl_powu_pre(y, le1, p, pi);
     864      154981 :     if (m != 1) break;
     865             :   }
     866      141243 :   *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      225463 : 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      225463 :   if (a==0) return a;
     879      225457 :   u2 = Fl_inv(l%r, r);
     880      225458 :   v = Fl_powu_pre(a, u2, p, pi);
     881      225456 :   w = Fl_powu_pre(v, l, p, pi);
     882      225456 :   w = pi? Fl_mul_pre(w, Fl_inv(a, p), p, pi): Fl_div(w, a, p);
     883      225446 :   if (w==1) return v;
     884      139170 :   if (y==0) y = Fl_lgener_pre_all(l, e, r, p, pi, &m);
     885      164584 :   while (w!=1)
     886             :   {
     887      144541 :     ulong k = 0;
     888      144541 :     p1 = w;
     889             :     do
     890             :     {
     891      188034 :       z = p1; p1 = Fl_powu_pre(p1, l, p, pi);
     892      188035 :       if (++k == e) return ULONG_MAX;
     893       68909 :     } while (p1!=1);
     894       25416 :     dl = Fl_log_pre(z, m, l, p, pi);
     895       25416 :     dl = Fl_neg(dl, l);
     896       25416 :     p1 = Fl_powu_pre(y,dl*upowuu(l,e-k-1),p,pi);
     897       25416 :     m = Fl_powu_pre(m, dl, p, pi);
     898       25416 :     e = k;
     899       25416 :     v = pi? Fl_mul_pre(p1,v,p,pi): Fl_mul(p1,v,p);
     900       25416 :     y = Fl_powu_pre(p1,l,p,pi);
     901       25416 :     w = pi? Fl_mul_pre(y,w,p,pi): Fl_mul(y,w,p);
     902             :   }
     903       20043 :   return v;
     904             : }
     905             : 
     906             : /* allow pi = 0 */
     907             : static ulong
     908      223104 : Fl_sqrtl_i(ulong a, ulong l, ulong p, ulong pi, ulong y, ulong m)
     909             : {
     910      223104 :   ulong r, e = u_lvalrem(p-1, l, &r);
     911      223104 :   return Fl_sqrtl_raw(a, l, e, r, p, pi, y, m);
     912             : }
     913             : /* allow pi = 0 */
     914             : ulong
     915      223104 : Fl_sqrtl_pre(ulong a, ulong l, ulong p, ulong pi)
     916      223104 : { 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      196109 : Fl_sqrtn_pre(ulong a, long n, ulong p, ulong pi, ulong *zetan)
     926             : {
     927      196109 :   ulong m, q = p-1, z;
     928      196109 :   ulong nn = n >= 0 ? (ulong)n: -(ulong)n;
     929      196109 :   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       79720 :   if (n==1)
     936             :   {
     937         420 :     if (zetan) *zetan = 1;
     938         420 :     return n < 0? Fl_inv(a,p): a;
     939             :   }
     940       79300 :   if (n==2)
     941             :   {
     942        6311 :     if (zetan) *zetan = p-1;
     943        6311 :     return Fl_sqrt_pre_i(a, 0, p, pi);
     944             :   }
     945       72989 :   if (a == 1 && !zetan) return a;
     946       43545 :   m = ugcd(nn, q);
     947       43545 :   z = 1;
     948       43545 :   if (m!=1)
     949             :   {
     950        2070 :     GEN F = factoru(m);
     951             :     long i, j, e;
     952             :     ulong r, zeta, y, l;
     953        4523 :     for (i = nbrows(F); i; i--)
     954             :     {
     955        2516 :       l = ucoeff(F,i,1);
     956        2516 :       j = ucoeff(F,i,2);
     957        2516 :       e = u_lvalrem(q,l, &r);
     958        2516 :       y = Fl_lgener_pre_all(l, e, r, p, pi, &zeta);
     959        2516 :       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        2516 :       if (a!=1)
     965             :         do
     966             :         {
     967        2355 :           a = Fl_sqrtl_raw(a, l, e, r, p, pi, y, zeta);
     968        2341 :           if (a==ULONG_MAX) return ULONG_MAX;
     969        2292 :         } while (--j);
     970             :     }
     971             :   }
     972       43482 :   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       43482 :   if (n < 0) a = Fl_inv(a, p);
     978       43482 :   if (zetan) *zetan = z;
     979       43482 :   return a;
     980             : }
     981             : 
     982             : ulong
     983      196109 : Fl_sqrtn(ulong a, long n, ulong p, ulong *zetan)
     984             : {
     985      196109 :   ulong pi = (p & HIGHMASK)? get_Fl_red(p): 0;
     986      196109 :   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       19198 : Fp_sqrts(long a, GEN p)
    1133             : {
    1134       19198 :   long v = vals(a)>>1;
    1135       19198 :   GEN r = gen_0;
    1136       19198 :   a >>= v << 1;
    1137       19198 :   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       15043 :     default:
    1185       15043 :       return gen_0;
    1186             :   }
    1187        2276 :   return remii(shifti(r, v), p);
    1188             : }
    1189             : 
    1190             : static GEN
    1191       77474 : Fp_sqrt_ii(GEN a, GEN y, GEN p)
    1192             : {
    1193       77474 :   pari_sp av = avma;
    1194       77474 :   GEN  q, v, w, p1 = subiu(p,1);
    1195       77477 :   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       77479 :   if (e == 1)
    1200             :   {
    1201       19091 :     q = addiu(shifti(p1,-2),1); /* (p+1) / 4 */
    1202       19088 :     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       19099 :     return e? v: NULL;
    1206             :   }
    1207       58388 :   as = itos_or_0(a);
    1208       58388 :   if (!as) as = itos_or_0(subii(a,p));
    1209       58392 :   if (as)
    1210             :   {
    1211       19198 :     GEN res = Fp_sqrts(as, p);
    1212       19198 :     if (!res) return gc_NULL(av);
    1213       19198 :     if (signe(res)) return gerepileupto(av, res);
    1214             :   }
    1215       56116 :   if (e == 2)
    1216             :   { /* Atkin's formula */
    1217       17284 :     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       17284 :     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       92605 :   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       53815 :     p1 = Fp_sqr(w,p);
    1244      110471 :     for (k=1; !equali1(p1) && k < e; k++) p1 = Fp_sqr(p1,p);
    1245       53816 :     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       53773 :     p1 = y;
    1248       76504 :     for (i=1; i < e-k; i++) p1 = Fp_sqr(p1,p);
    1249       53773 :     y = Fp_sqr(p1, p); e = k;
    1250       53773 :     w = Fp_mul(y, w, p);
    1251       53773 :     v = Fp_mul(v, p1, p);
    1252       53773 :     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     4890927 : Fp_sqrt_i(GEN a, GEN y, GEN p)
    1265             : {
    1266     4890927 :   pari_sp av = avma, av2;
    1267             :   GEN q;
    1268             : 
    1269     4890927 :   if (lgefint(p) == 3)
    1270             :   {
    1271     4813347 :     ulong pp = uel(p,2), u = umodiu(a, pp);
    1272     4813357 :     if (!u) return gen_0;
    1273     3602196 :     u = Fl_sqrt(u, pp);
    1274     3602300 :     return (u == ~0UL)? NULL: utoipos(u);
    1275             :   }
    1276       77580 :   a = modii(a, p); if (!signe(a)) return gen_0;
    1277       77474 :   a = Fp_sqrt_ii(a, y, p); if (!a) return gc_NULL(av);
    1278             :   /* smallest square root */
    1279       77077 :   av2 = avma; q = subii(p, a);
    1280       77074 :   if (cmpii(a, q) > 0) a = q; else set_avma(av2);
    1281       77076 :   return gerepileuptoint(av, a);
    1282             : }
    1283             : GEN
    1284     4833916 : Fp_sqrt(GEN a, GEN p) { return Fp_sqrt_i(a, NULL, p); }
    1285             : 
    1286             : /*********************************************************************/
    1287             : /**                        GCD & BEZOUT                             **/
    1288             : /*********************************************************************/
    1289             : 
    1290             : GEN
    1291    53600191 : lcmii(GEN x, GEN y)
    1292             : {
    1293             :   pari_sp av;
    1294             :   GEN a, b;
    1295    53600191 :   if (!signe(x) || !signe(y)) return gen_0;
    1296    53600202 :   av = avma; a = gcdii(x,y);
    1297    53599008 :   if (absequalii(a,y)) { set_avma(av); return absi(x); }
    1298    11945474 :   if (!equali1(a)) y = diviiexact(y,a);
    1299    11945466 :   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     5708710 : Fp_invgen(GEN x, GEN N, GEN *pd)
    1306             : {
    1307             :   GEN d, d0, e, v;
    1308     5708710 :   if (lgefint(N) == 3)
    1309             :   {
    1310     4965564 :     ulong dd, NN = N[2], xx = umodiu(x,NN);
    1311     4965612 :     if (!xx) { *pd = N; return gen_0; }
    1312     4965612 :     xx = Fl_invgen(xx, NN, &dd);
    1313     4966728 :     *pd = utoi(dd); return utoi(xx);
    1314             :   }
    1315      743146 :   *pd = d = bezout(x, N, &v, NULL);
    1316      743151 :   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      649358 :   e = diviiexact(N,d);
    1319      649358 :   d0 = Z_ppo(d, e); /* d = d0 d1, d0 coprime to N/d, rad(d1) | N/d */
    1320      649358 :   if (equali1(d0)) return v;
    1321      512229 :   if (!equalii(d,d0)) e = lcmii(e, diviiexact(d,d0));
    1322      512229 :   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     2414805 : gen_chinese(GEN x, GEN(*f)(GEN,GEN))
    1344             : {
    1345     2414805 :   GEN z = gassoc_proto(f,x,NULL);
    1346     2414798 :   if (z == gen_1) retmkintmod(gen_0,gen_1);
    1347     2414763 :   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      271300 : Z_chinese_pre(GEN A, GEN B, GEN *pC, GEN *pU, GEN *pd)
    1436             : {
    1437      271300 :   GEN u, d = bezout(A,B,&u,NULL); /* U = u(A/d), u(A/d) + v(B/d) = 1 */
    1438      271300 :   GEN t = diviiexact(A,d);
    1439      271288 :   *pU = mulii(u, t);
    1440      271279 :   *pC = mulii(t, B);
    1441      271284 :   if (pd) *pd = d;
    1442      271284 : }
    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     2976529 : Z_chinese_post(GEN a, GEN b, GEN C, GEN U, GEN d)
    1448             : {
    1449             :   GEN b_a;
    1450     2976529 :   if (!signe(a))
    1451             :   {
    1452      795904 :     if (d && !dvdii(b, d)) return NULL;
    1453      795904 :     return Fp_mul(b, U, C);
    1454             :   }
    1455     2180625 :   b_a = subii(b,a);
    1456     2180625 :   if (d && !dvdii(b_a, d)) return NULL;
    1457     2180625 :   return modii(addii(a, mulii(U, b_a)), C);
    1458             : }
    1459             : static ulong
    1460     1508583 : u_chinese_post(ulong a, ulong b, ulong C, ulong U)
    1461             : {
    1462     1508583 :   if (!a) return Fl_mul(b, U, C);
    1463     1508583 :   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      267596 : Z_chinese_all(GEN a, GEN b, GEN A, GEN B, GEN *pC)
    1475             : {
    1476      267596 :   GEN U; Z_chinese_pre(A, B, pC, &U, NULL);
    1477      267578 :   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      513489 : Z_chinese_coprime(GEN a, GEN b, GEN A, GEN B, GEN C)
    1484             : {
    1485      513489 :   pari_sp av = avma;
    1486      513489 :   GEN U = mulii(Fp_inv(A,B), A);
    1487      513489 :   return gerepileuptoint(av, Z_chinese_post(a,b,C,U, NULL));
    1488             : }
    1489             : ulong
    1490     1508471 : u_chinese_coprime(ulong a, ulong b, ulong A, ulong B, ulong C)
    1491     1508471 : { 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     2191493 : chinese1_coprime_Z_aux(GEN x, GEN y)
    1496             : {
    1497     2191493 :   GEN z = cgetg(3, t_INTMOD);
    1498     2191493 :   GEN A = gel(x,1), a = gel(x, 2);
    1499     2191493 :   GEN B = gel(y,1), b = gel(y, 2), C = mulii(A,B);
    1500     2191493 :   pari_sp av = avma;
    1501     2191493 :   GEN U = mulii(Fp_inv(A,B), A);
    1502     2191493 :   gel(z,2) = gerepileuptoint(av, Z_chinese_post(a,b,C,U, NULL));
    1503     2191493 :   gel(z,1) = C; return z;
    1504             : }
    1505             : GEN
    1506     2412460 : 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     2402907 : ZV_producttree(GEN xa)
    1514             : {
    1515     2402907 :   long n = lg(xa)-1;
    1516     2402907 :   long m = n==1 ? 1: expu(n-1)+1;
    1517     2402907 :   GEN T = cgetg(m+1, t_VEC), t;
    1518             :   long i, j, k;
    1519     2402904 :   t = cgetg(((n+1)>>1)+1, t_VEC);
    1520     2402904 :   if (typ(xa)==t_VECSMALL)
    1521             :   {
    1522     3133638 :     for (j=1, k=1; k<n; j++, k+=2)
    1523     2010550 :       gel(t, j) = muluu(xa[k], xa[k+1]);
    1524     1123088 :     if (k==n) gel(t, j) = utoi(xa[k]);
    1525             :   } else {
    1526     2649502 :     for (j=1, k=1; k<n; j++, k+=2)
    1527     1369683 :       gel(t, j) = mulii(gel(xa,k), gel(xa,k+1));
    1528     1279819 :     if (k==n) gel(t, j) = icopy(gel(xa,k));
    1529             :   }
    1530     2402906 :   gel(T,1) = t;
    1531     3797084 :   for (i=2; i<=m; i++)
    1532             :   {
    1533     1394171 :     GEN u = gel(T, i-1);
    1534     1394171 :     long n = lg(u)-1;
    1535     1394171 :     t = cgetg(((n+1)>>1)+1, t_VEC);
    1536     3132553 :     for (j=1, k=1; k<n; j++, k+=2)
    1537     1738375 :       gel(t, j) = mulii(gel(u, k), gel(u, k+1));
    1538     1394178 :     if (k==n) gel(t, j) = gel(u, k);
    1539     1394178 :     gel(T, i) = t;
    1540             :   }
    1541     2402913 :   return T;
    1542             : }
    1543             : 
    1544             : /* return [A mod P[i], i=1..#P], T = ZV_producttree(P) */
    1545             : GEN
    1546    57951166 : Z_ZV_mod_tree(GEN A, GEN P, GEN T)
    1547             : {
    1548             :   long i,j,k;
    1549    57951166 :   long m = lg(T)-1, n = lg(P)-1;
    1550             :   GEN t;
    1551    57951166 :   GEN Tp = cgetg(m+1, t_VEC);
    1552    57908970 :   gel(Tp, m) = mkvec(modii(A, gmael(T,m,1)));
    1553   122521203 :   for (i=m-1; i>=1; i--)
    1554             :   {
    1555    64692524 :     GEN u = gel(T, i);
    1556    64692524 :     GEN v = gel(Tp, i+1);
    1557    64692524 :     long n = lg(u)-1;
    1558    64692524 :     t = cgetg(n+1, t_VEC);
    1559   159344892 :     for (j=1, k=1; k<n; j++, k+=2)
    1560             :     {
    1561    94740425 :       gel(t, k)   = modii(gel(v, j), gel(u, k));
    1562    94769846 :       gel(t, k+1) = modii(gel(v, j), gel(u, k+1));
    1563             :     }
    1564    64604467 :     if (k==n) gel(t, k) = gel(v, j);
    1565    64604467 :     gel(Tp, i) = t;
    1566             :   }
    1567             :   {
    1568    57828679 :     GEN u = gel(T, i+1);
    1569    57828679 :     GEN v = gel(Tp, i+1);
    1570    57828679 :     long l = lg(u)-1;
    1571    57828679 :     if (typ(P)==t_VECSMALL)
    1572             :     {
    1573    55429972 :       GEN R = cgetg(n+1, t_VECSMALL);
    1574   205119460 :       for (j=1, k=1; j<=l; j++, k+=2)
    1575             :       {
    1576   149397950 :         uel(R,k) = umodiu(gel(v, j), P[k]);
    1577   149678936 :         if (k < n)
    1578   118491701 :           uel(R,k+1) = umodiu(gel(v, j), P[k+1]);
    1579             :       }
    1580    55721510 :       return R;
    1581             :     }
    1582             :     else
    1583             :     {
    1584     2398707 :       GEN R = cgetg(n+1, t_VEC);
    1585     6540519 :       for (j=1, k=1; j<=l; j++, k+=2)
    1586             :       {
    1587     4137959 :         gel(R,k) = modii(gel(v, j), gel(P,k));
    1588     4137961 :         if (k < n)
    1589     3377001 :           gel(R,k+1) = modii(gel(v, j), gel(P,k+1));
    1590             :       }
    1591     2402560 :       return R;
    1592             :     }
    1593             :   }
    1594             : }
    1595             : 
    1596             : /* T = ZV_producttree(P), R = ZV_chinesetree(P,T) */
    1597             : GEN
    1598    39366242 : ZV_chinese_tree(GEN A, GEN P, GEN T, GEN R)
    1599             : {
    1600    39366242 :   long m = lg(T)-1, n = lg(A)-1;
    1601             :   long i,j,k;
    1602    39366242 :   GEN Tp = cgetg(m+1, t_VEC);
    1603    39353603 :   GEN M = gel(T, 1);
    1604    39353603 :   GEN t = cgetg(lg(M), t_VEC);
    1605    39302833 :   if (typ(P)==t_VECSMALL)
    1606             :   {
    1607    88311591 :     for (j=1, k=1; k<n; j++, k+=2)
    1608             :     {
    1609    64913541 :       pari_sp av = avma;
    1610    64913541 :       GEN a = mului(A[k], gel(R,k)), b = mului(A[k+1], gel(R,k+1));
    1611    64805152 :       GEN tj = modii(addii(mului(P[k],b), mului(P[k+1],a)), gel(M,j));
    1612    64879729 :       gel(t, j) = gerepileuptoint(av, tj);
    1613             :     }
    1614    23398050 :     if (k==n) gel(t, j) = modii(mului(A[k], gel(R,k)), gel(M, j));
    1615             :   } else
    1616             :   {
    1617    33969675 :     for (j=1, k=1; k<n; j++, k+=2)
    1618             :     {
    1619    18020940 :       pari_sp av = avma;
    1620    18020940 :       GEN a = mulii(gel(A,k), gel(R,k)), b = mulii(gel(A,k+1), gel(R,k+1));
    1621    18023353 :       GEN tj = modii(addii(mulii(gel(P,k),b), mulii(gel(P,k+1),a)), gel(M,j));
    1622    18058714 :       gel(t, j) = gerepileuptoint(av, tj);
    1623             :     }
    1624    15948735 :     if (k==n) gel(t, j) = modii(mulii(gel(A,k), gel(R,k)), gel(M, j));
    1625             :   }
    1626    39340304 :   gel(Tp, 1) = t;
    1627    74725913 :   for (i=2; i<=m; i++)
    1628             :   {
    1629    35362611 :     GEN u = gel(T, i-1), M = gel(T, i);
    1630    35362611 :     GEN t = cgetg(lg(M), t_VEC);
    1631    35359440 :     GEN v = gel(Tp, i-1);
    1632    35359440 :     long n = lg(v)-1;
    1633    95025988 :     for (j=1, k=1; k<n; j++, k+=2)
    1634             :     {
    1635    59640379 :       pari_sp av = avma;
    1636    59589065 :       gel(t, j) = gerepileuptoint(av, modii(addii(mulii(gel(u, k), gel(v, k+1)),
    1637    59640379 :             mulii(gel(u, k+1), gel(v, k))), gel(M, j)));
    1638             :     }
    1639    35385609 :     if (k==n) gel(t, j) = gel(v, k);
    1640    35385609 :     gel(Tp, i) = t;
    1641             :   }
    1642    39363302 :   return gmael(Tp,m,1);
    1643             : }
    1644             : 
    1645             : static GEN
    1646     1441345 : ncV_polint_center_tree(GEN vA, GEN P, GEN T, GEN R, GEN m2)
    1647             : {
    1648     1441345 :   long i, l = lg(gel(vA,1)), n = lg(P);
    1649     1441345 :   GEN mod = gmael(T, lg(T)-1, 1), V = cgetg(l, t_COL);
    1650    34844227 :   for (i=1; i < l; i++)
    1651             :   {
    1652    33403163 :     pari_sp av = avma;
    1653    33403163 :     GEN c, A = cgetg(n, typ(P));
    1654             :     long j;
    1655   198815639 :     for (j=1; j < n; j++) A[j] = mael(vA,j,i);
    1656    33372082 :     c = Fp_center(ZV_chinese_tree(A, P, T, R), mod, m2);
    1657    33393870 :     gel(V,i) = gerepileuptoint(av, c);
    1658             :   }
    1659     1441064 :   return V;
    1660             : }
    1661             : 
    1662             : static GEN
    1663      638821 : nxV_polint_center_tree(GEN vA, GEN P, GEN T, GEN R, GEN m2)
    1664             : {
    1665      638821 :   long i, j, l, n = lg(P);
    1666      638821 :   GEN mod = gmael(T, lg(T)-1, 1), V, w;
    1667      638821 :   w = cgetg(n, t_VECSMALL);
    1668     2189837 :   for(j=1; j<n; j++) w[j] = lg(gel(vA,j));
    1669      638804 :   l = vecsmall_max(w);
    1670      638808 :   V = cgetg(l, t_POL);
    1671      638772 :   V[1] = mael(vA,1,1);
    1672     3903514 :   for (i=2; i < l; i++)
    1673             :   {
    1674     3264705 :     pari_sp av = avma;
    1675     3264705 :     GEN c, A = cgetg(n, typ(P));
    1676     3264284 :     if (typ(P)==t_VECSMALL)
    1677     7783640 :       for (j=1; j < n; j++) A[j] = i < w[j] ? mael(vA,j,i): 0;
    1678             :     else
    1679     4421702 :       for (j=1; j < n; j++) gel(A,j) = i < w[j] ? gmael(vA,j,i): gen_0;
    1680     3264284 :     c = Fp_center(ZV_chinese_tree(A, P, T, R), mod, m2);
    1681     3264761 :     gel(V,i) = gerepileuptoint(av, c);
    1682             :   }
    1683      638809 :   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       90863 :   for (i=1; i < l; i++)
    1693             :   {
    1694      334990 :     for (j=1; j < n; j++) gel(A,j) = gmael(vA,j,i);
    1695       86250 :     gel(V,i) = nxV_polint_center_tree(A, P, T, R, m2);
    1696             :   }
    1697        4613 :   return V;
    1698             : }
    1699             : 
    1700             : static GEN
    1701      365689 : polint_chinese(GEN worker, GEN mA, GEN P)
    1702             : {
    1703      365689 :   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      365689 :   if (l == 1) return cgetg(1, t_MAT);
    1709      336637 :   cnt = pending = 0;
    1710      336637 :   n = lg(P);
    1711      336637 :   A = cgetg(n, t_VEC);
    1712      336637 :   va = mkvec(A);
    1713      336637 :   M = cgetg(l, t_MAT);
    1714      336637 :   if (DEBUGLEVEL>4) timer_start(&ti);
    1715      336637 :   if (DEBUGLEVEL>5) err_printf("Start parallel Chinese remainder: ");
    1716      336637 :   mt_queue_start_lim(&pt, worker, l-1);
    1717     1300560 :   for (i=1; i<l || pending; i++)
    1718             :   {
    1719             :     long workid;
    1720     3661692 :     for(j=1; j < n; j++) gel(A,j) = gmael(mA,j,i);
    1721      963923 :     mt_queue_submit(&pt, i, i<l? va: NULL);
    1722      963923 :     done = mt_queue_get(&pt, &workid, &pending);
    1723      963923 :     if (done)
    1724             :     {
    1725      925053 :       gel(M,workid) = done;
    1726      925053 :       if (DEBUGLEVEL>5) err_printf("%ld%% ",(++cnt)*100/(l-1));
    1727             :     }
    1728             :   }
    1729      336637 :   if (DEBUGLEVEL>5) err_printf("\n");
    1730      336637 :   if (DEBUGLEVEL>4) timer_printf(&ti, "nmV_chinese_center");
    1731      336637 :   mt_queue_end(&pt);
    1732      336637 :   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      108605 : nmV_polint_center_tree_seq(GEN vA, GEN P, GEN T, GEN R, GEN m2)
    1764             : {
    1765      108605 :   long i, j, l = lg(gel(vA,1)), n = lg(P);
    1766      108605 :   GEN A = cgetg(n, t_VEC);
    1767      108605 :   GEN V = cgetg(l, t_MAT);
    1768      610039 :   for (i=1; i < l; i++)
    1769             :   {
    1770     2905136 :     for (j=1; j < n; j++) gel(A,j) = gmael(vA,j,i);
    1771      501435 :     gel(V,i) = ncV_polint_center_tree(A, P, T, R, m2);
    1772             :   }
    1773      108604 :   return V;
    1774             : }
    1775             : 
    1776             : GEN
    1777      924184 : nmV_polint_center_tree_worker(GEN vA, GEN T, GEN R, GEN P, GEN m2)
    1778             : {
    1779      924184 :   return ncV_polint_center_tree(vA, P, T, R, m2);
    1780             : }
    1781             : 
    1782             : static GEN
    1783      365598 : nmV_polint_center_tree(GEN mA, GEN P, GEN T, GEN R, GEN m2)
    1784             : {
    1785      365598 :   GEN worker = snm_closure(is_entry("_polint_worker"), mkvec4(T, R, P, m2));
    1786      365599 :   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     1423432 : ZX_nv_mod_tree(GEN B, GEN A, GEN T)
    1806             : {
    1807             :   pari_sp av;
    1808     1423432 :   long i, j, l = lg(B), n = lg(A)-1;
    1809     1423432 :   GEN V = cgetg(n+1, t_VEC);
    1810     6992667 :   for (j=1; j <= n; j++)
    1811             :   {
    1812     5569511 :     gel(V, j) = cgetg(l, t_VECSMALL);
    1813     5569251 :     mael(V, j, 1) = B[1]&VARNBITS;
    1814             :   }
    1815     1423156 :   av = avma;
    1816    12628455 :   for (i=2; i < l; i++)
    1817             :   {
    1818    11206386 :     GEN v = Z_ZV_mod_tree(gel(B, i), A, T);
    1819    78175463 :     for (j=1; j <= n; j++)
    1820    66978293 :       mael(V, j, i) = v[j];
    1821    11197170 :     set_avma(av);
    1822             :   }
    1823     6991857 :   for (j=1; j <= n; j++)
    1824     5569776 :     (void) Flx_renormalize(gel(V, j), l);
    1825     1422081 :   return V;
    1826             : }
    1827             : 
    1828             : static GEN
    1829      300382 : to_ZX(GEN a, long v) { return typ(a)==t_INT? scalarpol(a,v): a; }
    1830             : 
    1831             : GEN
    1832       31167 : ZXX_nv_mod_tree(GEN P, GEN xa, GEN T, long w)
    1833             : {
    1834       31167 :   pari_sp av = avma;
    1835       31167 :   long i, j, l = lg(P), n = lg(xa)-1;
    1836       31167 :   GEN V = cgetg(n+1, t_VEC);
    1837      126726 :   for (j=1; j <= n; j++)
    1838             :   {
    1839       95559 :     gel(V, j) = cgetg(l, t_POL);
    1840       95559 :     mael(V, j, 1) = P[1]&VARNBITS;
    1841             :   }
    1842      249825 :   for (i=2; i < l; i++)
    1843             :   {
    1844      218658 :     GEN v = ZX_nv_mod_tree(to_ZX(gel(P, i), w), xa, T);
    1845      869193 :     for (j=1; j <= n; j++)
    1846      650535 :       gmael(V, j, i) = gel(v,j);
    1847             :   }
    1848      126726 :   for (j=1; j <= n; j++)
    1849       95559 :     (void) FlxX_renormalize(gel(V, j), l);
    1850       31167 :   return gerepilecopy(av, V);
    1851             : }
    1852             : 
    1853             : GEN
    1854        4071 : ZXC_nv_mod_tree(GEN C, GEN xa, GEN T, long w)
    1855             : {
    1856        4071 :   pari_sp av = avma;
    1857        4071 :   long i, j, l = lg(C), n = lg(xa)-1;
    1858        4071 :   GEN V = cgetg(n+1, t_VEC);
    1859       17095 :   for (j = 1; j <= n; j++)
    1860       13024 :     gel(V, j) = cgetg(l, t_COL);
    1861       85794 :   for (i = 1; i < l; i++)
    1862             :   {
    1863       81722 :     GEN v = ZX_nv_mod_tree(to_ZX(gel(C, i), w), xa, T);
    1864      362572 :     for (j = 1; j <= n; j++)
    1865      280849 :       gmael(V, j, i) = gel(v,j);
    1866             :   }
    1867        4072 :   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     1195399 : ZV_nv_mod_tree(GEN B, GEN A, GEN T)
    1889             : {
    1890             :   pari_sp av;
    1891     1195399 :   long i, j, l = lg(B), n = lg(A)-1;
    1892     1195399 :   GEN V = cgetg(n+1, t_VEC);
    1893     6103736 :   for (j=1; j <= n; j++) gel(V, j) = cgetg(l, t_VECSMALL);
    1894     1195330 :   av = avma;
    1895    45525537 :   for (i=1; i < l; i++)
    1896             :   {
    1897    44334331 :     GEN v = Z_ZV_mod_tree(gel(B, i), A, T);
    1898   246596986 :     for (j=1; j <= n; j++) mael(V, j, i) = v[j];
    1899    44297112 :     set_avma(av);
    1900             :   }
    1901     1191206 :   return V;
    1902             : }
    1903             : 
    1904             : GEN
    1905      175315 : ZM_nv_mod_tree(GEN M, GEN xa, GEN T)
    1906             : {
    1907      175315 :   pari_sp av = avma;
    1908      175315 :   long i, j, l = lg(M), n = lg(xa)-1;
    1909      175315 :   GEN V = cgetg(n+1, t_VEC);
    1910      863066 :   for (j=1; j <= n; j++) gel(V, j) = cgetg(l, t_MAT);
    1911     1370460 :   for (i=1; i < l; i++)
    1912             :   {
    1913     1195147 :     GEN v = ZV_nv_mod_tree(gel(M, i), xa, T);
    1914     6102866 :     for (j=1; j <= n; j++) gmael(V, j, i) = gel(v,j);
    1915             :   }
    1916      175313 :   return gerepilecopy(av, V);
    1917             : }
    1918             : 
    1919             : static GEN
    1920     2398700 : ZV_sqr(GEN z)
    1921             : {
    1922     2398700 :   long i,l = lg(z);
    1923     2398700 :   GEN x = cgetg(l, t_VEC);
    1924     2398698 :   if (typ(z)==t_VECSMALL)
    1925     5543745 :     for (i=1; i<l; i++) gel(x,i) = sqru(z[i]);
    1926             :   else
    1927     4344319 :     for (i=1; i<l; i++) gel(x,i) = sqri(gel(z,i));
    1928     2398677 :   return x;
    1929             : }
    1930             : 
    1931             : static GEN
    1932    12266267 : ZT_sqr(GEN x)
    1933             : {
    1934    12266267 :   if (typ(x) == t_INT) return sqri(x);
    1935    16052717 :   pari_APPLY_type(t_VEC, ZT_sqr(gel(x,i)))
    1936             : }
    1937             : 
    1938             : static GEN
    1939     2398695 : ZV_invdivexact(GEN y, GEN x)
    1940             : {
    1941     2398695 :   long i, l = lg(y);
    1942     2398695 :   GEN z = cgetg(l,t_VEC);
    1943     2398695 :   if (typ(x)==t_VECSMALL)
    1944     5543478 :     for (i=1; i<l; i++)
    1945             :     {
    1946     4420751 :       pari_sp av = avma;
    1947     4420751 :       ulong a = Fl_inv(umodiu(diviuexact(gel(y,i),x[i]), x[i]), x[i]);
    1948     4420949 :       set_avma(av); gel(z,i) = utoi(a);
    1949             :     }
    1950             :   else
    1951     4344323 :     for (i=1; i<l; i++)
    1952     3068356 :       gel(z,i) = Fp_inv(diviiexact(gel(y,i), gel(x,i)), gel(x,i));
    1953     2398694 :   return z;
    1954             : }
    1955             : 
    1956             : /* P t_VECSMALL or t_VEC of t_INT  */
    1957             : GEN
    1958     2398692 : ZV_chinesetree(GEN P, GEN T)
    1959             : {
    1960     2398692 :   GEN T2 = ZT_sqr(T), P2 = ZV_sqr(P);
    1961     2398695 :   GEN mod = gmael(T,lg(T)-1,1);
    1962     2398695 :   return ZV_invdivexact(Z_ZV_mod_tree(mod, P2, T2), P);
    1963             : }
    1964             : 
    1965             : static GEN
    1966      926660 : gc_chinese(pari_sp av, GEN T, GEN a, GEN *pt_mod)
    1967             : {
    1968      926660 :   if (!pt_mod)
    1969       12377 :     return gerepileupto(av, a);
    1970             :   else
    1971             :   {
    1972      914283 :     GEN mod = gmael(T, lg(T)-1, 1);
    1973      914283 :     gerepileall(av, 2, &a, &mod);
    1974      914285 :     *pt_mod = mod;
    1975      914285 :     return a;
    1976             :   }
    1977             : }
    1978             : 
    1979             : GEN
    1980      152844 : ZV_chinese_center(GEN A, GEN P, GEN *pt_mod)
    1981             : {
    1982      152844 :   pari_sp av = avma;
    1983      152844 :   GEN T = ZV_producttree(P);
    1984      152844 :   GEN R = ZV_chinesetree(P, T);
    1985      152844 :   GEN a = ZV_chinese_tree(A, P, T, R);
    1986      152844 :   GEN mod = gmael(T, lg(T)-1, 1);
    1987      152844 :   GEN ca = Fp_center(a, mod, shifti(mod,-1));
    1988      152844 :   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      159802 : nxV_chinese_center_tree(GEN A, GEN P, GEN T, GEN R)
    2003             : {
    2004      159802 :   pari_sp av = avma;
    2005      159802 :   GEN m2 = shifti(gmael(T, lg(T)-1, 1), -1);
    2006      159802 :   GEN a = nxV_polint_center_tree(A, P, T, R, m2);
    2007      159803 :   return gerepileupto(av, a);
    2008             : }
    2009             : 
    2010             : GEN
    2011      392756 : nxV_chinese_center(GEN A, GEN P, GEN *pt_mod)
    2012             : {
    2013      392756 :   pari_sp av = avma;
    2014      392756 :   GEN T = ZV_producttree(P);
    2015      392756 :   GEN R = ZV_chinesetree(P, T);
    2016      392754 :   GEN m2 = shifti(gmael(T, lg(T)-1, 1), -1);
    2017      392756 :   GEN a = nxV_polint_center_tree(A, P, T, R, m2);
    2018      392756 :   return gc_chinese(av, T, a, pt_mod);
    2019             : }
    2020             : 
    2021             : GEN
    2022       10263 : ncV_chinese_center(GEN A, GEN P, GEN *pt_mod)
    2023             : {
    2024       10263 :   pari_sp av = avma;
    2025       10263 :   GEN T = ZV_producttree(P);
    2026       10263 :   GEN R = ZV_chinesetree(P, T);
    2027       10263 :   GEN m2 = shifti(gmael(T, lg(T)-1, 1), -1);
    2028       10263 :   GEN a = ncV_polint_center_tree(A, P, T, R, m2);
    2029       10263 :   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      108605 : nmV_chinese_center_tree_seq(GEN A, GEN P, GEN T, GEN R)
    2052             : {
    2053      108605 :   pari_sp av = avma;
    2054      108605 :   GEN m2 = shifti(gmael(T, lg(T)-1, 1), -1);
    2055      108605 :   GEN a = nmV_polint_center_tree_seq(A, P, T, R, m2);
    2056      108604 :   return gerepileupto(av, a);
    2057             : }
    2058             : 
    2059             : GEN
    2060      365599 : nmV_chinese_center(GEN A, GEN P, GEN *pt_mod)
    2061             : {
    2062      365599 :   pari_sp av = avma;
    2063      365599 :   GEN T = ZV_producttree(P);
    2064      365599 :   GEN R = ZV_chinesetree(P, T);
    2065      365598 :   GEN m2 = shifti(gmael(T, lg(T)-1, 1), -1);
    2066      365598 :   GEN a = nmV_polint_center_tree(A, P, T, R, m2);
    2067      365599 :   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    11791272 : Fl_2powu_pre(ulong n, ulong p, ulong pi)
    2117             : {
    2118    11791272 :   ulong y = 2;
    2119    11791272 :   int j = 1+bfffo(n);
    2120             :   /* normalize, i.e set highest bit to 1 (we know n != 0) */
    2121    11791272 :   n<<=j; j = BITS_IN_LONG-j; /* first bit is now implicit */
    2122   516672193 :   for (; j; n<<=1,j--)
    2123             :   {
    2124   504910060 :     y = Fl_sqr_pre(y,p,pi);
    2125   504900942 :     if (n & HIGHBIT) y = Fl_double(y, p);
    2126             :   }
    2127    11762133 :   return y;
    2128             : }
    2129             : 
    2130             : /* 2^n mod p; assume n > 1 and !(p & HIGHMASK) */
    2131             : static ulong
    2132     4289389 : Fl_2powu(ulong n, ulong p)
    2133             : {
    2134     4289389 :   ulong y = 2;
    2135     4289389 :   int j = 1+bfffo(n);
    2136             :   /* normalize, i.e set highest bit to 1 (we know n != 0) */
    2137     4289389 :   n<<=j; j = BITS_IN_LONG-j; /* first bit is now implicit */
    2138    25361414 :   for (; j; n<<=1,j--)
    2139             :   {
    2140    21071979 :     y = (y*y) % p;
    2141    21071979 :     if (n & HIGHBIT) y = Fl_double(y, p);
    2142             :   }
    2143     4289435 :   return y;
    2144             : }
    2145             : 
    2146             : /* allow pi = 0 */
    2147             : ulong
    2148    97967370 : Fl_powu_pre(ulong x, ulong n0, ulong p, ulong pi)
    2149             : {
    2150             :   ulong y, z, n;
    2151    97967370 :   if (!pi) return Fl_powu(x, n0, p);
    2152    95527193 :   if (n0 <= 1)
    2153             :   { /* frequent special cases */
    2154     6806931 :     if (n0 == 1) return x;
    2155      104245 :     if (n0 == 0) return 1;
    2156             :   }
    2157    88720280 :   if (x <= 2)
    2158             :   {
    2159    12048492 :     if (x == 2) return Fl_2powu_pre(n0, p, pi);
    2160      256898 :     return x; /* 0 or 1 */
    2161             :   }
    2162    76671788 :   y = 1; z = x; n = n0;
    2163             :   for(;;)
    2164             :   {
    2165   543959453 :     if (n&1) y = Fl_mul_pre(y,z,p,pi);
    2166   543822234 :     n>>=1; if (!n) return y;
    2167   467017901 :     z = Fl_sqr_pre(z,p,pi);
    2168             :   }
    2169             : }
    2170             : 
    2171             : ulong
    2172   140141915 : Fl_powu(ulong x, ulong n0, ulong p)
    2173             : {
    2174             :   ulong y, z, n;
    2175   140141915 :   if (n0 <= 2)
    2176             :   { /* frequent special cases */
    2177    65713125 :     if (n0 == 2) return Fl_sqr(x,p);
    2178    32127824 :     if (n0 == 1) return x;
    2179     1912293 :     if (n0 == 0) return 1;
    2180             :   }
    2181    74398836 :   if (x <= 1) return x; /* 0 or 1 */
    2182    73965603 :   if (p & HIGHMASK)
    2183     7880999 :     return Fl_powu_pre(x, n0, p, get_Fl_red(p));
    2184    66084604 :   if (x == 2) return Fl_2powu(n0, p);
    2185    61795214 :   y = 1; z = x; n = n0;
    2186             :   for(;;)
    2187             :   {
    2188   274952912 :     if (n&1) y = (y*z) % p;
    2189   274952912 :     n>>=1; if (!n) return y;
    2190   213157698 :     z = (z*z) % p;
    2191             :   }
    2192             : }
    2193             : 
    2194             : /* Reduce data dependency to maximize internal parallelism; allow pi = 0 */
    2195             : GEN
    2196    12788190 : Fl_powers_pre(ulong x, long n, ulong p, ulong pi)
    2197             : {
    2198             :   long i, k;
    2199    12788190 :   GEN z = cgetg(n + 2, t_VECSMALL);
    2200    12778635 :   z[1] = 1; if (n == 0) return z;
    2201    12778635 :   z[2] = x;
    2202    12778635 :   if (pi)
    2203             :   {
    2204    90041371 :     for (i = 3, k=2; i <= n; i+=2, k++)
    2205             :     {
    2206    77467575 :       z[i] = Fl_sqr_pre(z[k], p, pi);
    2207    77471757 :       z[i+1] = Fl_mul_pre(z[k], z[k+1], p, pi);
    2208             :     }
    2209    12573796 :     if (i==n+1) z[i] = Fl_sqr_pre(z[k], p, pi);
    2210             :   }
    2211      212878 :   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   400522881 :     for (i = 2; i <= n; i++) z[i+1] = (z[i] * x) % p;
    2222    12788987 :   return z;
    2223             : }
    2224             : 
    2225             : GEN
    2226      295924 : Fl_powers(ulong x, long n, ulong p)
    2227             : {
    2228      295924 :   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       14994 : Fp_invmBarrett(GEN p, long s)
    2245             : {
    2246       14994 :   GEN R, Q = dvmdii(int2n(3*s),p,&R);
    2247       14994 :   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     9212246 : Fp_rem_mBarrett(GEN a, GEN B, long s, GEN N)
    2254             : {
    2255     9212246 :   pari_sp av = avma;
    2256     9212246 :   GEN P = gel(B, 1), Q = gel(B, 2); /* 2^(3s) = P N + Q, 0 <= Q < N */
    2257     9212246 :   long t = expi(P)+1; /* 2^(t-1) <= P < 2^t */
    2258     9212246 :   GEN u = shifti(a, -3*s), v = remi2n(a, 3*s); /* a = 2^(3s)u + v */
    2259     9212246 :   GEN A = addii(v, mulii(Q,u)); /* 0 <= A < 2^(3s+1) */
    2260     9212246 :   GEN q = shifti(mulii(shifti(A, t-3*s), P), -t); /* A/N - 4 < q <= A/N */
    2261     9212246 :   GEN r = subii(A, mulii(q, N));
    2262     9212246 :   GEN sr= subii(r,N);     /* 0 <= r < 4*N */
    2263     9212246 :   if (signe(sr)<0) return gerepileuptoint(av, r);
    2264     5028279 :   r=sr; sr = subii(r,N);  /* 0 <= r < 3*N */
    2265     5028279 :   if (signe(sr)<0) return gerepileuptoint(av, r);
    2266      186280 :   r=sr; sr = subii(r,N);  /* 0 <= r < 2*N */
    2267      186280 :   return gerepileuptoint(av, signe(sr)>=0 ? sr:r);
    2268             : }
    2269             : 
    2270             : /* Montgomery reduction */
    2271             : 
    2272             : INLINE ulong
    2273      697920 : 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    69007718 : _sqr_montred(void * E, GEN x)
    2284             : {
    2285    69007718 :   struct montred * D = (struct montred *) E;
    2286    69007718 :   return red_montgomery(sqri(x), D->N, D->inv);
    2287             : }
    2288             : 
    2289             : /* Montgomery reduction */
    2290             : static GEN
    2291     7070372 : _mul_montred(void * E, GEN x, GEN y)
    2292             : {
    2293     7070372 :   struct montred * D = (struct montred *) E;
    2294     7070372 :   return red_montgomery(mulii(x, y), D->N, D->inv);
    2295             : }
    2296             : 
    2297             : static GEN
    2298    11252599 : _mul2_montred(void * E, GEN x)
    2299             : {
    2300    11252599 :   struct montred * D = (struct montred *) E;
    2301    11252599 :   GEN z = shifti(_sqr_montred(E, x), 1);
    2302    11250887 :   long l = lgefint(D->N);
    2303    11888711 :   while (lgefint(z) > l) z = subii(z, D->N);
    2304    11251145 :   return z;
    2305             : }
    2306             : 
    2307             : static GEN
    2308    23203742 : _sqr_remii(void* N, GEN x)
    2309    23203742 : { return remii(sqri(x), (GEN) N); }
    2310             : 
    2311             : static GEN
    2312     1520741 : _mul_remii(void* N, GEN x, GEN y)
    2313     1520741 : { return remii(mulii(x, y), (GEN) N); }
    2314             : 
    2315             : static GEN
    2316     3228594 : _mul2_remii(void* N, GEN x)
    2317     3228594 : { 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     8437940 : _sqr_remiibar(void *E, GEN x)
    2327             : {
    2328     8437940 :   struct redbarrett * D = (struct redbarrett *) E;
    2329     8437940 :   return Fp_rem_mBarrett(sqri(x), D->iM, D->s, D->N);
    2330             : }
    2331             : 
    2332             : static GEN
    2333      774306 : _mul_remiibar(void *E, GEN x, GEN y)
    2334             : {
    2335      774306 :   struct redbarrett * D = (struct redbarrett *) E;
    2336      774306 :   return Fp_rem_mBarrett(mulii(x, y), D->iM, D->s, D->N);
    2337             : }
    2338             : 
    2339             : static GEN
    2340     2079690 : _mul2_remiibar(void *E, GEN x)
    2341             : {
    2342     2079690 :   struct redbarrett * D = (struct redbarrett *) E;
    2343     2079690 :   return Fp_double(_sqr_remiibar(E, x), D->N);
    2344             : }
    2345             : 
    2346             : static long
    2347      893674 : Fp_select_red(GEN *y, ulong k, GEN N, long lN, muldata *D, void **pt_E)
    2348             : {
    2349      893674 :   if (lN >= Fp_POW_BARRETT_LIMIT && (k==0 || ((double)k)*expi(*y) > 2 + expi(N)))
    2350             :   {
    2351       14994 :     struct redbarrett * E = (struct redbarrett *) stack_malloc(sizeof(struct redbarrett));
    2352       14994 :     D->sqr = &_sqr_remiibar;
    2353       14994 :     D->mul = &_mul_remiibar;
    2354       14994 :     D->mul2 = &_mul2_remiibar;
    2355       14994 :     E->N = N;
    2356       14994 :     E->s = 1+(expi(N)>>1);
    2357       14994 :     E->iM = Fp_invmBarrett(N, E->s);
    2358       14994 :     *pt_E = (void*) E;
    2359       14994 :     return 0;
    2360             :   }
    2361      878680 :   else if (mod2(N) && lN < Fp_POW_REDC_LIMIT)
    2362             :   {
    2363      697916 :     struct montred * E = (struct montred *) stack_malloc(sizeof(struct montred));
    2364      697913 :     *y = remii(shifti(*y, bit_accuracy(lN)), N);
    2365      697921 :     D->sqr = &_sqr_montred;
    2366      697921 :     D->mul = &_mul_montred;
    2367      697921 :     D->mul2 = &_mul2_montred;
    2368      697921 :     E->N = N;
    2369      697921 :     E->inv = init_montdata(N);
    2370      697917 :     *pt_E = (void*) E;
    2371      697917 :     return 1;
    2372             :   }
    2373             :   else
    2374             :   {
    2375      180773 :     D->sqr = &_sqr_remii;
    2376      180773 :     D->mul = &_mul_remii;
    2377      180773 :     D->mul2 = &_mul2_remii;
    2378      180773 :     *pt_E = (void*) N;
    2379      180773 :     return 0;
    2380             :   }
    2381             : }
    2382             : 
    2383             : GEN
    2384     3017315 : Fp_powu(GEN A, ulong k, GEN N)
    2385             : {
    2386     3017315 :   long lN = lgefint(N);
    2387             :   int base_is_2, use_montgomery;
    2388             :   muldata D;
    2389             :   void *E;
    2390             :   pari_sp av;
    2391             : 
    2392     3017315 :   if (lN == 3) {
    2393     1540080 :     ulong n = uel(N,2);
    2394     1540080 :     return utoi( Fl_powu(umodiu(A, n), k, n) );
    2395             :   }
    2396     1477235 :   if (k <= 2)
    2397             :   { /* frequent special cases */
    2398      854468 :     if (k == 2) return Fp_sqr(A,N);
    2399      300220 :     if (k == 1) return A;
    2400           0 :     if (k == 0) return gen_1;
    2401             :   }
    2402      622767 :   av = avma; A = modii(A,N);
    2403      622767 :   base_is_2 = 0;
    2404      622767 :   if (lgefint(A) == 3) switch(A[2])
    2405             :   {
    2406         861 :     case 1: set_avma(av); return gen_1;
    2407       43832 :     case 2:  base_is_2 = 1; break;
    2408             :   }
    2409             : 
    2410             :   /* TODO: Move this out of here and use for general modular computations */
    2411      621906 :   use_montgomery = Fp_select_red(&A, k, N, lN, &D, &E);
    2412      621907 :   if (base_is_2)
    2413       43832 :     A = gen_powu_fold_i(A, k, E, D.sqr, D.mul2);
    2414             :   else
    2415      578075 :     A = gen_powu_i(A, k, E, D.sqr, D.mul);
    2416      621907 :   if (use_montgomery)
    2417             :   {
    2418      527248 :     A = red_montgomery(A, N, ((struct montred *) E)->inv);
    2419      527248 :     if (cmpii(A, N) >= 0) A = subii(A,N);
    2420             :   }
    2421      621907 :   return gerepileuptoint(av, A);
    2422             : }
    2423             : 
    2424             : GEN
    2425       29383 : Fp_pows(GEN A, long k, GEN N)
    2426             : {
    2427       29383 :   if (lgefint(N) == 3) {
    2428       14925 :     ulong n = N[2];
    2429       14925 :     ulong a = umodiu(A, n);
    2430       14925 :     if (k < 0) {
    2431         133 :       a = Fl_inv(a, n);
    2432         133 :       k = -k;
    2433             :     }
    2434       14925 :     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    36224868 : Fp_pow(GEN A, GEN K, GEN N)
    2443             : {
    2444             :   pari_sp av;
    2445    36224868 :   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    36224868 :   s = signe(K);
    2452    36224868 :   if (!s) return dvdii(A,N)? gen_0: gen_1;
    2453    35194096 :   if (lN == 3 && lgefint(K) == 3)
    2454             :   {
    2455    34448238 :     ulong n = N[2], a = umodiu(A, n);
    2456    34448521 :     if (s < 0) a = Fl_inv(a, n);
    2457    34448603 :     if (a <= 1) return utoi(a); /* 0 or 1 */
    2458    30904505 :     return utoi(Fl_powu(a, uel(K,2), n));
    2459             :   }
    2460             : 
    2461      745858 :   av = avma;
    2462      745858 :   if (s < 0) y = Fp_inv(A,N);
    2463             :   else
    2464             :   {
    2465      743921 :     y = modii(A,N);
    2466      744112 :     if (!signe(y)) { set_avma(av); return gen_0; }
    2467             :   }
    2468      746049 :   if (lgefint(K) == 3) return gerepileuptoint(av, Fp_powu(y, K[2], N));
    2469             : 
    2470      271975 :   base_is_2 = 0;
    2471      271975 :   sy = abscmpii(y, shifti(N,-1)) > 0;
    2472      271964 :   if (sy) y = subii(N,y);
    2473      271970 :   sA = sy && mod2(K);
    2474      271969 :   if (lgefint(y) == 3) switch(y[2])
    2475             :   {
    2476         207 :     case 1:  set_avma(av); return sA ? subis(N,1): gen_1;
    2477      147017 :     case 2:  base_is_2 = 1; break;
    2478             :   }
    2479             : 
    2480             :   /* TODO: Move this out of here and use for general modular computations */
    2481      271762 :   use_montgomery = Fp_select_red(&y, 0UL, N, lN, &D, &E);
    2482      271784 :   if (base_is_2)
    2483      147032 :     y = gen_pow_fold_i(y, K, E, D.sqr, D.mul2);
    2484             :   else
    2485      124752 :     y = gen_pow_i(y, K, E, D.sqr, D.mul);
    2486      271810 :   if (use_montgomery)
    2487             :   {
    2488      170685 :     y = red_montgomery(y, N, ((struct montred *) E)->inv);
    2489      170685 :     if (cmpii(y,N) >= 0) y = subii(y,N);
    2490             :   }
    2491      271808 :   if (sA) y = subii(N, y);
    2492      271807 :   return gerepileuptoint(av,y);
    2493             : }
    2494             : 
    2495             : static GEN
    2496    14129270 : _Fp_mul(void *E, GEN x, GEN y) { return Fp_mul(x,y,(GEN)E); }
    2497             : static GEN
    2498     8134253 : _Fp_sqr(void *E, GEN x) { return Fp_sqr(x,(GEN)E); }
    2499             : static GEN
    2500       47162 : _Fp_one(void *E) { (void) E; return gen_1; }
    2501             : 
    2502             : GEN
    2503         105 : Fp_pow_init(GEN x, GEN n, long k, GEN p)
    2504         105 : { return gen_pow_init(x, n, k, (void*)p, &_Fp_sqr, &_Fp_mul); }
    2505             : 
    2506             : GEN
    2507       43694 : Fp_pow_table(GEN R, GEN n, GEN p)
    2508       43694 : { return gen_pow_table(R, n, (void*)p, &_Fp_one, &_Fp_mul); }
    2509             : 
    2510             : GEN
    2511        5931 : Fp_powers(GEN x, long n, GEN p)
    2512             : {
    2513        5931 :   if (lgefint(p) == 3)
    2514        2463 :     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    27875651 : _Fp_pow(void *E, GEN x, GEN n) { return Fp_pow(x,n,(GEN)E); }
    2523             : static GEN
    2524         147 : _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      927370 : _Fp_red(void *E, GEN x) { return Fp_red(x, (GEN)E); }
    2532             : static GEN
    2533     1221959 : _Fp_add(void *E, GEN x, GEN y) { (void) E; return addii(x,y); }
    2534             : static GEN
    2535     1123230 : _Fp_neg(void *E, GEN x) { (void) E; return negi(x); }
    2536             : static GEN
    2537      609429 : _Fp_rmul(void *E, GEN x, GEN y) { (void) E; return mulii(x,y); }
    2538             : static GEN
    2539       36647 : _Fp_inv(void *E, GEN x) { return Fp_inv(x,(GEN)E); }
    2540             : static int
    2541      265848 : _Fp_equal0(GEN x) { return signe(x)==0; }
    2542             : static GEN
    2543       33797 : _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        7870 : const struct bb_field *get_Fp_field(void **E, GEN p)
    2549        7870 : { *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       61476 : Fp_log_halfgcd(ulong bnd, GEN C, GEN g, GEN p)
    2648             : {
    2649       61476 :   pari_sp av = avma;
    2650             :   GEN h1, h2, F, G;
    2651       61476 :   if (!Fp_ratlift(g,p,C,shifti(C,-1),&h1,&h2)) return gc_NULL(av);
    2652       37036 :   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       36910 :   return gc_NULL(av);
    2660             : }
    2661             : 
    2662             : static GEN
    2663       61476 : Fp_log_find_rel(GEN b, ulong bnd, GEN C, GEN p, GEN *g, long *e)
    2664             : {
    2665             :   GEN rel;
    2666       61476 :   do { (*e)++; *g = Fp_mul(*g, b, p); rel = Fp_log_halfgcd(bnd, C, *g, p); }
    2667       61476 :   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             : static long
    2679       59731 : tr(long i) { return odd(i) ? (i+1)>>1: -(i>>1); }
    2680             : 
    2681             : static long
    2682      169813 : rt(long i) { return i>0 ? 2*i-1: -2*i; }
    2683             : 
    2684             : /* add u^e */
    2685             : static void
    2686        2163 : addifsmooth1(struct Fp_log_rel *r, GEN z, long u, long e)
    2687             : {
    2688        2163 :   pari_sp av = avma;
    2689        2163 :   long off = r->prmax+1;
    2690        2163 :   GEN F = cgetg(3, t_MAT);
    2691        2163 :   gel(F,1) = vecsmall_append(gel(z,1), off+rt(u));
    2692        2163 :   gel(F,2) = vecsmall_append(gel(z,2), e);
    2693        2163 :   gel(r->rel,++r->nbrel) = gerepileupto(av, F);
    2694        2163 : }
    2695             : 
    2696             : /* add u^-1 v^-1 */
    2697             : static void
    2698       83825 : addifsmooth2(struct Fp_log_rel *r, GEN z, long u, long v)
    2699             : {
    2700       83825 :   pari_sp av = avma;
    2701       83825 :   long off = r->prmax+1;
    2702       83825 :   GEN P = mkvecsmall2(off+rt(u),off+rt(v)), E = mkvecsmall2(-1,-1);
    2703       83825 :   GEN F = cgetg(3, t_MAT);
    2704       83825 :   gel(F,1) = vecsmall_concat(gel(z,1), P);
    2705       83825 :   gel(F,2) = vecsmall_concat(gel(z,2), E);
    2706       83825 :   gel(r->rel,++r->nbrel) = gerepileupto(av, F);
    2707       83825 : }
    2708             : 
    2709             : /*
    2710             : Let p=C^2+c
    2711             : Solve h = (C+x)*(C+a)-p = 0 [mod l]
    2712             : h= -c+x*(C+a)+C*a = 0  [mod l]
    2713             : x = (c-C*a)/(C+a) [mod l]
    2714             : h = -c+C*(x+a)+a*x
    2715             : */
    2716             : 
    2717             : GEN
    2718       30251 : Fp_log_sieve_worker(long a, long prmax, GEN C, GEN c, GEN Ci, GEN ci, GEN pi, GEN sz)
    2719             : {
    2720       30251 :   pari_sp ltop = avma;
    2721       30251 :   long i, j, th, n = lg(pi)-1, rel = 1, ab = labs(a), ae;
    2722       30251 :   GEN sieve = zero_zv(2*ab+2)+1+ab;
    2723       30257 :   GEN L = cgetg(1+2*ab+2, t_VEC);
    2724       30252 :   pari_sp av = avma;
    2725       30252 :   GEN z, h = addis(C,a);
    2726       30240 :   if ((z = Z_issmooth_fact(h, prmax)))
    2727             :   {
    2728        2169 :     gel(L, rel++) = mkvec2(z, mkvecsmall3(1, a, -1));
    2729        2168 :     av = avma;
    2730             :   }
    2731    12461654 :   for (i=1; i<=n; i++)
    2732             :   {
    2733    12432604 :     ulong li = pi[i], s = sz[i], al = smodss(a,li);
    2734    12416181 :     ulong iv = Fl_invsafe(Fl_add(Ci[i],al,li),li);
    2735             :     long u;
    2736    12687654 :     if (!iv) continue;
    2737    12373462 :     u = Fl_add(Fl_mul(Fl_sub(ci[i],Fl_mul(Ci[i],al,li),li), iv ,li),ab%li,li)-ab;
    2738    46183849 :     for(j = u; j<=ab; j+=li) sieve[j] += s;
    2739             :   }
    2740       29050 :   if (a)
    2741             :   {
    2742       30148 :     long e = expi(mulis(C,a));
    2743       30180 :     th = e - expu(e) - 1;
    2744          54 :   } else th = -1;
    2745       30254 :   ae = a>=0 ? ab-1: ab;
    2746    15518004 :   for (j = 1-ab; j <= ae; j++)
    2747    15487431 :     if (sieve[j]>=th)
    2748             :     {
    2749      108843 :       GEN h = absi(addis(subii(mulis(C,a+j),c), a*j));
    2750      108712 :       if ((z = Z_issmooth_fact(h, prmax)))
    2751             :       {
    2752      106448 :         gel(L, rel++) = mkvec2(z, mkvecsmall3(2, a, j));
    2753      106478 :         av = avma;
    2754        2289 :       } else set_avma(av);
    2755             :     }
    2756             :   /* j = a */
    2757       30573 :   if (sieve[a]>=th)
    2758             :   {
    2759         448 :     GEN h = absi(addiu(subii(mulis(C,2*a),c), a*a));
    2760         448 :     if ((z = Z_issmooth_fact(h, prmax)))
    2761         364 :       gel(L, rel++) = mkvec2(z, mkvecsmall3(1, a, -2));
    2762             :   }
    2763       30573 :   setlg(L, rel); return gerepilecopy(ltop, L);
    2764             : }
    2765             : 
    2766             : static long
    2767          63 : Fp_log_sieve(struct Fp_log_rel *r, GEN C, GEN c, GEN Ci, GEN ci, GEN pi, GEN sz)
    2768             : {
    2769             :   struct pari_mt pt;
    2770          63 :   long i, j, nb = 0;
    2771          63 :   GEN worker = snm_closure(is_entry("_Fp_log_sieve_worker"),
    2772             :                mkvecn(7, utoi(r->prmax), C, c, Ci, ci, pi, sz));
    2773          63 :   long running, pending = 0;
    2774          63 :   GEN W = zerovec(r->nbgen);
    2775          63 :   mt_queue_start_lim(&pt, worker, r->nbgen);
    2776       30459 :   for (i = 0; (running = (i < r->nbgen)) || pending; i++)
    2777             :   {
    2778             :     GEN done;
    2779             :     long idx;
    2780       30396 :     mt_queue_submit(&pt, i, running ? mkvec(stoi(tr(i))): NULL);
    2781       30396 :     done = mt_queue_get(&pt, &idx, &pending);
    2782       30396 :     if (!done || lg(done)==1) continue;
    2783       27636 :     gel(W, idx+1) = done;
    2784       27636 :     nb += lg(done)-1;
    2785       27636 :     if (DEBUGLEVEL && (i&127)==0)
    2786           0 :       err_printf("%ld%% ",100*nb/r->nbmax);
    2787             :   }
    2788          63 :   mt_queue_end(&pt);
    2789       26362 :   for(j = 1; j <= r->nbgen && r->nbrel < r->nbmax; j++)
    2790             :   {
    2791             :     long ll, m;
    2792       26299 :     GEN L = gel(W,j);
    2793       26299 :     if (isintzero(L)) continue;
    2794       23681 :     ll = lg(L);
    2795      109669 :     for (m=1; m<ll && r->nbrel < r->nbmax ; m++)
    2796             :     {
    2797       85988 :       GEN Lm = gel(L,m), h = gel(Lm, 1), v = gel(Lm, 2);
    2798       85988 :       if (v[1] == 1)
    2799        2163 :         addifsmooth1(r, h, v[2], v[3]);
    2800             :       else
    2801       83825 :         addifsmooth2(r, h, v[2], v[3]);
    2802             :     }
    2803             :   }
    2804          63 :   return j;
    2805             : }
    2806             : 
    2807             : static GEN
    2808         837 : ECP_psi(GEN x, GEN y)
    2809             : {
    2810         837 :   long prec = realprec(x);
    2811         837 :   GEN lx = glog(x, prec), ly = glog(y, prec);
    2812         837 :   GEN u = gdiv(lx, ly);
    2813         837 :   return gpow(u, gneg(u),prec);
    2814             : }
    2815             : 
    2816             : struct computeG
    2817             : {
    2818             :   GEN C;
    2819             :   long bnd, nbi;
    2820             : };
    2821             : 
    2822             : static GEN
    2823         837 : _computeG(void *E, GEN gen)
    2824             : {
    2825         837 :   struct computeG * d = (struct computeG *) E;
    2826         837 :   GEN ps = ECP_psi(gmul(gen,d->C), stoi(d->bnd));
    2827         837 :   return gsub(gmul(gsqr(gen),ps),gmulgs(gaddgs(gen,d->nbi),3));
    2828             : }
    2829             : 
    2830             : static long
    2831          63 : compute_nbgen(GEN C, long bnd, long nbi)
    2832             : {
    2833             :   struct computeG d;
    2834          63 :   d.C = shifti(C, 1);
    2835          63 :   d.bnd = bnd;
    2836          63 :   d.nbi = nbi;
    2837          63 :   return itos(ground(zbrent((void*)&d, _computeG, gen_2, stoi(bnd), DEFAULTPREC)));
    2838             : }
    2839             : 
    2840             : static GEN
    2841        1714 : _psi(void*E, GEN y)
    2842             : {
    2843        1714 :   GEN lx = (GEN) E;
    2844        1714 :   long prec = realprec(lx);
    2845        1714 :   GEN ly = glog(y, prec);
    2846        1714 :   GEN u = gdiv(lx, ly);
    2847        1714 :   return gsub(gdiv(y ,ly), gpow(u, u, prec));
    2848             : }
    2849             : 
    2850             : static GEN
    2851          63 : opt_param(GEN x, long prec)
    2852             : {
    2853          63 :   return zbrent((void*)glog(x,prec), _psi, gen_2, x, prec);
    2854             : }
    2855             : 
    2856             : static GEN
    2857          63 : check_kernel(long nbg, long N, long prmax, GEN C, GEN M, GEN p, GEN m)
    2858             : {
    2859          63 :   pari_sp av = avma;
    2860          63 :   long lM = lg(M)-1, nbcol = lM;
    2861          63 :   long tbs = maxss(1, expu(nbg/expi(m)));
    2862             :   for (;;)
    2863          42 :   {
    2864         105 :     GEN K = FpMs_leftkernel_elt_col(M, nbcol, N, m);
    2865             :     GEN tab;
    2866         105 :     long i, f=0;
    2867         105 :     long l = lg(K), lm = lgefint(m);
    2868         105 :     GEN idx = diviiexact(subiu(p,1),m), g;
    2869             :     pari_timer ti;
    2870         105 :     if (DEBUGLEVEL) timer_start(&ti);
    2871         210 :     for(i=1; i<l; i++)
    2872         210 :       if (signe(gel(K,i)))
    2873         105 :         break;
    2874         105 :     g = Fp_pow(utoi(i), idx, p);
    2875         105 :     tab = Fp_pow_init(g, p, tbs, p);
    2876         105 :     K = FpC_Fp_mul(K, Fp_inv(gel(K,i), m), m);
    2877      121520 :     for(i=1; i<l; i++)
    2878             :     {
    2879      121415 :       GEN k = gel(K,i);
    2880      121415 :       GEN j = i<=prmax ? utoi(i): addis(C,tr(i-(prmax+1)));
    2881      121415 :       if (signe(k)==0 || !equalii(Fp_pow_table(tab, k, p), Fp_pow(j, idx, p)))
    2882       82369 :         gel(K,i) = cgetineg(lm);
    2883             :       else
    2884       39046 :         f++;
    2885             :     }
    2886         105 :     if (DEBUGLEVEL) timer_printf(&ti,"found %ld/%ld logs", f, nbg);
    2887         105 :     if(f > (nbg>>1)) return gerepileupto(av, K);
    2888       10024 :     for(i=1; i<=nbcol; i++)
    2889             :     {
    2890        9982 :       long a = 1+random_Fl(lM);
    2891        9982 :       swap(gel(M,a),gel(M,i));
    2892             :     }
    2893          42 :     if (4*nbcol>5*nbg) nbcol = nbcol*9/10;
    2894             :   }
    2895             : }
    2896             : 
    2897             : static GEN
    2898         126 : Fp_log_find_ind(GEN a, GEN K, long prmax, GEN C, GEN p, GEN m)
    2899             : {
    2900         126 :   pari_sp av=avma;
    2901         126 :   GEN aa = gen_1;
    2902         126 :   long AV = 0;
    2903             :   for(;;)
    2904           0 :   {
    2905         126 :     GEN A = Fp_log_find_rel(a, prmax, C, p, &aa, &AV);
    2906         126 :     GEN F = gel(A,1), E = gel(A,2);
    2907         126 :     GEN Ao = gen_0;
    2908         126 :     long i, l = lg(F);
    2909         771 :     for(i=1; i<l; i++)
    2910             :     {
    2911         645 :       GEN Ki = gel(K,F[i]);
    2912         645 :       if (signe(Ki)<0) break;
    2913         645 :       Ao = addii(Ao, mulis(Ki, E[i]));
    2914             :     }
    2915         126 :     if (i==l) return Fp_divu(Ao, AV, m);
    2916           0 :     aa = gerepileuptoint(av, aa);
    2917             :   }
    2918             : }
    2919             : 
    2920             : static GEN
    2921          63 : Fp_log_index(GEN a, GEN b, GEN m, GEN p)
    2922             : {
    2923          63 :   pari_sp av = avma, av2;
    2924          63 :   long i, j, nbi, nbr = 0, nbrow, nbg;
    2925             :   GEN C, c, Ci, ci, pi, pr, sz, l, Ao, Bo, K, d, p_1;
    2926             :   pari_timer ti;
    2927             :   struct Fp_log_rel r;
    2928          63 :   ulong bnds = itou(roundr_safe(opt_param(sqrti(p),DEFAULTPREC)));
    2929          63 :   ulong bnd = 4*bnds;
    2930          63 :   if (!bnds || cmpii(sqru(bnds),m)>=0) return NULL;
    2931             : 
    2932          63 :   p_1 = subiu(p,1);
    2933          63 :   if (!is_pm1(gcdii(m,diviiexact(p_1,m))))
    2934           0 :     m = diviiexact(p_1, Z_ppo(p_1, m));
    2935          63 :   pr = primes_upto_zv(bnd);
    2936          63 :   nbi = lg(pr)-1;
    2937          63 :   C = sqrtremi(p, &c);
    2938          63 :   av2 = avma;
    2939       12796 :   for (i = 1; i <= nbi; ++i)
    2940             :   {
    2941       12733 :     ulong lp = pr[i];
    2942       26894 :     while (lp <= bnd)
    2943             :     {
    2944       14161 :       nbr++;
    2945       14161 :       lp *= pr[i];
    2946             :     }
    2947             :   }
    2948          63 :   pi = cgetg(nbr+1,t_VECSMALL);
    2949          63 :   Ci = cgetg(nbr+1,t_VECSMALL);
    2950          63 :   ci = cgetg(nbr+1,t_VECSMALL);
    2951          63 :   sz = cgetg(nbr+1,t_VECSMALL);
    2952       12796 :   for (i = 1, j = 1; i <= nbi; ++i)
    2953             :   {
    2954       12733 :     ulong lp = pr[i], sp = expu(2*lp-1);
    2955       26894 :     while (lp <= bnd)
    2956             :     {
    2957       14161 :       pi[j] = lp;
    2958       14161 :       Ci[j] = umodiu(C, lp);
    2959       14161 :       ci[j] = umodiu(c, lp);
    2960       14161 :       sz[j] = sp;
    2961       14161 :       lp *= pr[i];
    2962       14161 :       j++;
    2963             :     }
    2964             :   }
    2965          63 :   r.nbrel = 0;
    2966          63 :   r.nbgen = compute_nbgen(C, bnd, nbi);
    2967          63 :   r.nbmax = 2*(nbi+r.nbgen);
    2968          63 :   r.rel = cgetg(r.nbmax+1,t_VEC);
    2969          63 :   r.prmax = pr[nbi];
    2970          63 :   if (DEBUGLEVEL)
    2971             :   {
    2972           0 :     err_printf("bnd=%lu Size FB=%ld extra gen=%ld \n", bnd, nbi, r.nbgen);
    2973           0 :     timer_start(&ti);
    2974             :   }
    2975          63 :   nbg = Fp_log_sieve(&r, C, c, Ci, ci, pi, sz);
    2976          63 :   nbrow = r.prmax + nbg;
    2977          63 :   if (DEBUGLEVEL)
    2978             :   {
    2979           0 :     err_printf("\n");
    2980           0 :     timer_printf(&ti," %ld relations, %ld generators", r.nbrel, nbi+nbg);
    2981             :   }
    2982          63 :   setlg(r.rel,r.nbrel+1);
    2983          63 :   r.rel = gerepilecopy(av2, r.rel);
    2984          63 :   K = check_kernel(nbi+nbrow-r.prmax, nbrow, r.prmax, C, r.rel, p, m);
    2985          63 :   if (DEBUGLEVEL) timer_start(&ti);
    2986          63 :   Ao = Fp_log_find_ind(a, K, r.prmax, C, p, m);
    2987          63 :   if (DEBUGLEVEL) timer_printf(&ti," log element");
    2988          63 :   Bo = Fp_log_find_ind(b, K, r.prmax, C, p, m);
    2989          63 :   if (DEBUGLEVEL) timer_printf(&ti," log generator");
    2990          63 :   d = gcdii(Ao,Bo);
    2991          63 :   l = Fp_div(diviiexact(Ao, d) ,diviiexact(Bo, d), m);
    2992          63 :   if (!equalii(a,Fp_pow(b,l,p))) pari_err_BUG("Fp_log_index");
    2993          63 :   return gerepileuptoint(av, l);
    2994             : }
    2995             : 
    2996             : static int
    2997     4663466 : Fp_log_use_index(long e, long p)
    2998             : {
    2999     4663466 :   return (e >= 27 && 20*(p+6)<=e*e);
    3000             : }
    3001             : 
    3002             : /* Trivial cases a = 1, -1. Return x s.t. g^x = a or [] if no such x exist */
    3003             : static GEN
    3004     8458579 : Fp_easylog(void *E, GEN a, GEN g, GEN ord)
    3005             : {
    3006     8458579 :   pari_sp av = avma;
    3007     8458579 :   GEN p = (GEN)E;
    3008             :   /* assume a reduced mod p, p not necessarily prime */
    3009     8458579 :   if (equali1(a)) return gen_0;
    3010             :   /* p > 2 */
    3011     5438740 :   if (equalii(subiu(p,1), a))  /* -1 */
    3012             :   {
    3013             :     pari_sp av2;
    3014             :     GEN t;
    3015     1321829 :     ord = get_arith_Z(ord);
    3016     1321829 :     if (mpodd(ord)) { set_avma(av); return cgetg(1, t_VEC); } /* no solution */
    3017     1321815 :     t = shifti(ord,-1); /* only possible solution */
    3018     1321815 :     av2 = avma;
    3019     1321815 :     if (!equalii(Fp_pow(g, t, p), a)) { set_avma(av); return cgetg(1, t_VEC); }
    3020     1321787 :     set_avma(av2); return gerepileuptoint(av, t);
    3021             :   }
    3022     4116912 :   if (typ(ord)==t_INT && BPSW_psp(p) && Fp_log_use_index(expi(ord),expi(p)))
    3023          63 :     return Fp_log_index(a, g, ord, p);
    3024     4116849 :   return gc_NULL(av); /* not easy */
    3025             : }
    3026             : 
    3027             : GEN
    3028     3925704 : Fp_log(GEN a, GEN g, GEN ord, GEN p)
    3029             : {
    3030     3925704 :   GEN v = get_arith_ZZM(ord);
    3031     3925668 :   GEN F = gmael(v,2,1);
    3032     3925668 :   long lF = lg(F)-1, lmax;
    3033     3925668 :   if (lF == 0) return equali1(a)? gen_0: cgetg(1, t_VEC);
    3034     3925640 :   lmax = expi(gel(F,lF));
    3035     3925639 :   if (BPSW_psp(p) && Fp_log_use_index(lmax,expi(p)))
    3036          91 :     v = mkvec2(gel(v,1),ZM_famat_limit(gel(v,2),int2n(27)));
    3037     3925631 :   return gen_PH_log(a,g,v,(void*)p,&Fp_star);
    3038             : }
    3039             : 
    3040             : /* assume !(p & HIGHMASK) */
    3041             : static ulong
    3042      131788 : Fl_log_naive(ulong a, ulong g, ulong ord, ulong p)
    3043             : {
    3044      131788 :   ulong i, h=1;
    3045      363417 :   for (i = 0; i < ord; i++, h = (h * g) % p)
    3046      363417 :     if (a==h) return i;
    3047           0 :   return ~0UL;
    3048             : }
    3049             : 
    3050             : static ulong
    3051       25017 : Fl_log_naive_pre(ulong a, ulong g, ulong ord, ulong p, ulong pi)
    3052             : {
    3053       25017 :   ulong i, h=1;
    3054       64141 :   for (i = 0; i < ord; i++, h = Fl_mul_pre(h, g, p, pi))
    3055       64141 :     if (a==h) return i;
    3056           0 :   return ~0UL;
    3057             : }
    3058             : 
    3059             : static ulong
    3060           0 : Fl_log_Fp(ulong a, ulong g, ulong ord, ulong p)
    3061             : {
    3062           0 :   pari_sp av = avma;
    3063           0 :   GEN r = Fp_log(utoi(a),utoi(g),utoi(ord),utoi(p));
    3064           0 :   return gc_ulong(av, typ(r)==t_INT ? itou(r): ~0UL);
    3065             : }
    3066             : 
    3067             : /* allow pi = 0 */
    3068             : ulong
    3069       25416 : Fl_log_pre(ulong a, ulong g, ulong ord, ulong p, ulong pi)
    3070             : {
    3071       25416 :   if (!pi) return Fl_log(a, g, ord, p);
    3072       25017 :   if (ord <= 200) return Fl_log_naive_pre(a, g, ord, p, pi);
    3073           0 :   return Fl_log_Fp(a, g, ord, p);
    3074             : }
    3075             : 
    3076             : ulong
    3077      131788 : Fl_log(ulong a, ulong g, ulong ord, ulong p)
    3078             : {
    3079      131788 :   if (ord <= 200)
    3080           0 :     return (p&HIGHMASK)? Fl_log_naive_pre(a, g, ord, p, get_Fl_red(p))
    3081      131788 :                        : Fl_log_naive(a, g, ord, p);
    3082           0 :   return Fl_log_Fp(a, g, ord, p);
    3083             : }
    3084             : 
    3085             : /* find x such that h = g^x mod N > 1, N = prod_{i <= l} P[i]^E[i], P[i] prime.
    3086             :  * PHI[l] = eulerphi(N / P[l]^E[l]).   Destroys P/E */
    3087             : static GEN
    3088         126 : znlog_rec(GEN h, GEN g, GEN N, GEN P, GEN E, GEN PHI)
    3089             : {
    3090         126 :   long l = lg(P) - 1, e = E[l];
    3091         126 :   GEN p = gel(P, l), phi = gel(PHI,l), pe = e == 1? p: powiu(p, e);
    3092             :   GEN a,b, hp,gp, hpe,gpe, ogpe; /* = order(g mod p^e) | p^(e-1)(p-1) */
    3093             : 
    3094         126 :   if (l == 1) {
    3095          98 :     hpe = h;
    3096          98 :     gpe = g;
    3097             :   } else {
    3098          28 :     hpe = modii(h, pe);
    3099          28 :     gpe = modii(g, pe);
    3100             :   }
    3101         126 :   if (e == 1) {
    3102          42 :     hp = hpe;
    3103          42 :     gp = gpe;
    3104             :   } else {
    3105          84 :     hp = remii(hpe, p);
    3106          84 :     gp = remii(gpe, p);
    3107             :   }
    3108         126 :   if (hp == gen_0 || gp == gen_0) return NULL;
    3109         105 :   if (absequaliu(p, 2))
    3110             :   {
    3111          35 :     GEN n = int2n(e);
    3112          35 :     ogpe = Zp_order(gpe, gen_2, e, n);
    3113          35 :     a = Fp_log(hpe, gpe, ogpe, n);
    3114          35 :     if (typ(a) != t_INT) return NULL;
    3115             :   }
    3116             :   else
    3117             :   { /* Avoid black box groups: (Z/p^2)^* / (Z/p)^* ~ (Z/pZ, +), where DL
    3118             :        is trivial */
    3119             :     /* [order(gp), factor(order(gp))] */
    3120          70 :     GEN v = Fp_factored_order(gp, subiu(p,1), p);
    3121          70 :     GEN ogp = gel(v,1);
    3122          70 :     if (!equali1(Fp_pow(hp, ogp, p))) return NULL;
    3123          70 :     a = Fp_log(hp, gp, v, p);
    3124          70 :     if (typ(a) != t_INT) return NULL;
    3125          70 :     if (e == 1) ogpe = ogp;
    3126             :     else
    3127             :     { /* find a s.t. g^a = h (mod p^e), p odd prime, e > 0, (h,p) = 1 */
    3128             :       /* use p-adic log: O(log p + e) mul*/
    3129             :       long vpogpe, vpohpe;
    3130             : 
    3131          28 :       hpe = Fp_mul(hpe, Fp_pow(gpe, negi(a), pe), pe);
    3132          28 :       gpe = Fp_pow(gpe, ogp, pe);
    3133             :       /* g,h = 1 mod p; compute b s.t. h = g^b */
    3134             : 
    3135             :       /* v_p(order g mod pe) */
    3136          28 :       vpogpe = equali1(gpe)? 0: e - Z_pval(subiu(gpe,1), p);
    3137             :       /* v_p(order h mod pe) */
    3138          28 :       vpohpe = equali1(hpe)? 0: e - Z_pval(subiu(hpe,1), p);
    3139          28 :       if (vpohpe > vpogpe) return NULL;
    3140             : 
    3141          28 :       ogpe = mulii(ogp, powiu(p, vpogpe)); /* order g mod p^e */
    3142          28 :       if (is_pm1(gpe)) return is_pm1(hpe)? a: NULL;
    3143          28 :       b = gdiv(Qp_log(cvtop(hpe, p, e)), Qp_log(cvtop(gpe, p, e)));
    3144          28 :       a = addii(a, mulii(ogp, padic_to_Q(b)));
    3145             :     }
    3146             :   }
    3147             :   /* gp^a = hp => x = a mod ogpe => generalized Pohlig-Hellman strategy */
    3148          91 :   if (l == 1) return a;
    3149             : 
    3150          28 :   N = diviiexact(N, pe); /* make N coprime to p */
    3151          28 :   h = Fp_mul(h, Fp_pow(g, modii(negi(a), phi), N), N);
    3152          28 :   g = Fp_pow(g, modii(ogpe, phi), N);
    3153          28 :   setlg(P, l); /* remove last element */
    3154          28 :   setlg(E, l);
    3155          28 :   b = znlog_rec(h, g, N, P, E, PHI);
    3156          28 :   if (!b) return NULL;
    3157          28 :   return addmulii(a, b, ogpe);
    3158             : }
    3159             : 
    3160             : static GEN
    3161          98 : get_PHI(GEN P, GEN E)
    3162             : {
    3163          98 :   long i, l = lg(P);
    3164          98 :   GEN PHI = cgetg(l, t_VEC);
    3165          98 :   gel(PHI,1) = gen_1;
    3166         126 :   for (i=1; i<l-1; i++)
    3167             :   {
    3168          28 :     GEN t, p = gel(P,i);
    3169          28 :     long e = E[i];
    3170          28 :     t = mulii(powiu(p, e-1), subiu(p,1));
    3171          28 :     if (i > 1) t = mulii(t, gel(PHI,i));
    3172          28 :     gel(PHI,i+1) = t;
    3173             :   }
    3174          98 :   return PHI;
    3175             : }
    3176             : 
    3177             : GEN
    3178         238 : znlog(GEN h, GEN g, GEN o)
    3179             : {
    3180         238 :   pari_sp av = avma;
    3181             :   GEN N, fa, P, E, x;
    3182         238 :   switch (typ(g))
    3183             :   {
    3184          28 :     case t_PADIC:
    3185             :     {
    3186          28 :       GEN p = gel(g,2);
    3187          28 :       long v = valp(g);
    3188          28 :       if (v < 0) pari_err_DIM("znlog");
    3189          28 :       if (v > 0) {
    3190           0 :         long k = gvaluation(h, p);
    3191           0 :         if (k % v) return cgetg(1,t_VEC);
    3192           0 :         k /= v;
    3193           0 :         if (!gequal(h, gpowgs(g,k))) { set_avma(av); return cgetg(1,t_VEC); }
    3194           0 :         return gc_stoi(av, k);
    3195             :       }
    3196          28 :       N = gel(g,3);
    3197          28 :       g = Rg_to_Fp(g, N);
    3198          28 :       break;
    3199             :     }
    3200         203 :     case t_INTMOD:
    3201         203 :       N = gel(g,1);
    3202         203 :       g = gel(g,2); break;
    3203           7 :     default: pari_err_TYPE("znlog", g);
    3204             :       return NULL; /* LCOV_EXCL_LINE */
    3205             :   }
    3206         231 :   if (equali1(N)) { set_avma(av); return gen_0; }
    3207         231 :   h = Rg_to_Fp(h, N);
    3208         224 :   if (o) return gerepileupto(av, Fp_log(h, g, o, N));
    3209          98 :   fa = Z_factor(N);
    3210          98 :   P = gel(fa,1);
    3211          98 :   E = vec_to_vecsmall(gel(fa,2));
    3212          98 :   x = znlog_rec(h, g, N, P, E, get_PHI(P,E));
    3213          98 :   if (!x) { set_avma(av); return cgetg(1,t_VEC); }
    3214          63 :   return gerepileuptoint(av, x);
    3215             : }
    3216             : 
    3217             : GEN
    3218      173306 : Fp_sqrtn(GEN a, GEN n, GEN p, GEN *zeta)
    3219             : {
    3220      173306 :   if (lgefint(p)==3)
    3221             :   {
    3222      172872 :     long nn = itos_or_0(n);
    3223      172872 :     if (nn)
    3224             :     {
    3225      172872 :       ulong pp = p[2];
    3226             :       ulong uz;
    3227      172872 :       ulong r = Fl_sqrtn(umodiu(a,pp),nn,pp, zeta ? &uz:NULL);
    3228      172851 :       if (r==ULONG_MAX) return NULL;
    3229      172795 :       if (zeta) *zeta = utoi(uz);
    3230      172795 :       return utoi(r);
    3231             :     }
    3232             :   }
    3233         434 :   a = modii(a,p);
    3234         434 :   if (!signe(a))
    3235             :   {
    3236           0 :     if (zeta) *zeta = gen_1;
    3237           0 :     if (signe(n) < 0) pari_err_INV("Fp_sqrtn", mkintmod(gen_0,p));
    3238           0 :     return gen_0;
    3239             :   }
    3240         434 :   if (absequaliu(n,2))
    3241             :   {
    3242         238 :     if (zeta) *zeta = subiu(p,1);
    3243         238 :     return signe(n) > 0 ? Fp_sqrt(a,p): Fp_sqrt(Fp_inv(a, p),p);
    3244             :   }
    3245         196 :   return gen_Shanks_sqrtn(a,n,subiu(p,1),zeta,(void*)p,&Fp_star);
    3246             : }
    3247             : 
    3248             : /*********************************************************************/
    3249             : /**                              FACTORIAL                          **/
    3250             : /*********************************************************************/
    3251             : GEN
    3252       89675 : mulu_interval_step(ulong a, ulong b, ulong step)
    3253             : {
    3254       89675 :   pari_sp av = avma;
    3255             :   ulong k, l, N, n;
    3256             :   long lx;
    3257             :   GEN x;
    3258             : 
    3259       89675 :   if (!a) return gen_0;
    3260       89675 :   if (step == 1) return mulu_interval(a, b);
    3261       89675 :   n = 1 + (b-a) / step;
    3262       89675 :   b -= (b-a) % step;
    3263       89675 :   if (n < 61)
    3264             :   {
    3265       88292 :     if (n == 1) return utoipos(a);
    3266       67911 :     x = muluu(a,a+step); if (n == 2) return x;
    3267      533952 :     for (k=a+2*step; k<=b; k+=step) x = mului(k,x);
    3268       53220 :     return gerepileuptoint(av, x);
    3269             :   }
    3270             :   /* step | b-a */
    3271        1383 :   lx = 1; x = cgetg(2 + n/2, t_VEC);
    3272        1384 :   N = b + a;
    3273        1384 :   for (k = a;; k += step)
    3274             :   {
    3275      227447 :     l = N - k; if (l <= k) break;
    3276      226063 :     gel(x,lx++) = muluu(k,l);
    3277             :   }
    3278        1384 :   if (l == k) gel(x,lx++) = utoipos(k);
    3279        1384 :   setlg(x, lx);
    3280        1384 :   return gerepileuptoint(av, ZV_prod(x));
    3281             : }
    3282             : /* return a * (a+1) * ... * b. Assume a <= b  [ note: factoring out powers of 2
    3283             :  * first is slower ... ] */
    3284             : GEN
    3285      158924 : mulu_interval(ulong a, ulong b)
    3286             : {
    3287      158924 :   pari_sp av = avma;
    3288             :   ulong k, l, N, n;
    3289             :   long lx;
    3290             :   GEN x;
    3291             : 
    3292      158924 :   if (!a) return gen_0;
    3293      158924 :   n = b - a + 1;
    3294      158924 :   if (n < 61)
    3295             :   {
    3296      158195 :     if (n == 1) return utoipos(a);
    3297      107879 :     x = muluu(a,a+1); if (n == 2) return x;
    3298      403873 :     for (k=a+2; k<b; k++) x = mului(k,x);
    3299             :     /* avoid k <= b: broken if b = ULONG_MAX */
    3300       93777 :     return gerepileuptoint(av, mului(b,x));
    3301             :   }
    3302         729 :   lx = 1; x = cgetg(2 + n/2, t_VEC);
    3303         733 :   N = b + a;
    3304         733 :   for (k = a;; k++)
    3305             :   {
    3306       27473 :     l = N - k; if (l <= k) break;
    3307       26743 :     gel(x,lx++) = muluu(k,l);
    3308             :   }
    3309         730 :   if (l == k) gel(x,lx++) = utoipos(k);
    3310         730 :   setlg(x, lx);
    3311         730 :   return gerepileuptoint(av, ZV_prod(x));
    3312             : }
    3313             : GEN
    3314         560 : muls_interval(long a, long b)
    3315             : {
    3316         560 :   pari_sp av = avma;
    3317         560 :   long lx, k, l, N, n = b - a + 1;
    3318             :   GEN x;
    3319             : 
    3320         560 :   if (a <= 0 && b >= 0) return gen_0;
    3321         287 :   if (n < 61)
    3322             :   {
    3323         287 :     x = stoi(a);
    3324         511 :     for (k=a+1; k<=b; k++) x = mulsi(k,x);
    3325         287 :     return gerepileuptoint(av, x);
    3326             :   }
    3327           0 :   lx = 1; x = cgetg(2 + n/2, t_VEC);
    3328           0 :   N = b + a;
    3329           0 :   for (k = a;; k++)
    3330             :   {
    3331           0 :     l = N - k; if (l <= k) break;
    3332           0 :     gel(x,lx++) = mulss(k,l);
    3333             :   }
    3334           0 :   if (l == k) gel(x,lx++) = stoi(k);
    3335           0 :   setlg(x, lx);
    3336           0 :   return gerepileuptoint(av, ZV_prod(x));
    3337             : }
    3338             : 
    3339             : GEN
    3340         105 : mpprimorial(long n)
    3341             : {
    3342         105 :   pari_sp av = avma;
    3343         105 :   if (n <= 12) switch(n)
    3344             :   {
    3345          14 :     case 0: case 1: return gen_1;
    3346           7 :     case 2: return gen_2;
    3347          14 :     case 3: case 4: return utoipos(6);
    3348          14 :     case 5: case 6: return utoipos(30);
    3349          28 :     case 7: case 8: case 9: case 10: return utoipos(210);
    3350          14 :     case 11: case 12: return utoipos(2310);
    3351           7 :     default: pari_err_DOMAIN("primorial", "argument","<",gen_0,stoi(n));
    3352             :   }
    3353           7 :   return gerepileuptoint(av, zv_prod_Z(primes_upto_zv(n)));
    3354             : }
    3355             : 
    3356             : GEN
    3357      496338 : mpfact(long n)
    3358             : {
    3359      496338 :   pari_sp av = avma;
    3360             :   GEN a, v;
    3361             :   long k;
    3362      496338 :   if (n <= 12) switch(n)
    3363             :   {
    3364      428633 :     case 0: case 1: return gen_1;
    3365       24331 :     case 2: return gen_2;
    3366        3388 :     case 3: return utoipos(6);
    3367        4145 :     case 4: return utoipos(24);
    3368        2887 :     case 5: return utoipos(120);
    3369        2542 :     case 6: return utoipos(720);
    3370        2448 :     case 7: return utoipos(5040);
    3371        2437 :     case 8: return utoipos(40320);
    3372        2458 :     case 9: return utoipos(362880);
    3373        2694 :     case 10:return utoipos(3628800);
    3374        1409 :     case 11:return utoipos(39916800);
    3375         577 :     case 12:return utoipos(479001600);
    3376           0 :     default: pari_err_DOMAIN("factorial", "argument","<",gen_0,stoi(n));
    3377             :   }
    3378       18390 :   v = cgetg(expu(n) + 2, t_VEC);
    3379       18381 :   for (k = 1;; k++)
    3380       85881 :   {
    3381      104262 :     long m = n >> (k-1), l;
    3382      104262 :     if (m <= 2) break;
    3383       85872 :     l = (1 + (n >> k)) | 1;
    3384             :     /* product of odd numbers in ]n / 2^k, 2 / 2^(k-1)] */
    3385       85872 :     a = mulu_interval_step(l, m, 2);
    3386       85830 :     gel(v,k) = k == 1? a: powiu(a, k);
    3387             :   }
    3388       85931 :   a = gel(v,--k); while (--k) a = mulii(a, gel(v,k));
    3389       18390 :   a = shifti(a, factorial_lval(n, 2));
    3390       18388 :   return gerepileuptoint(av, a);
    3391             : }
    3392             : 
    3393             : ulong
    3394       56679 : factorial_Fl(long n, ulong p)
    3395             : {
    3396             :   long k;
    3397             :   ulong v;
    3398       56679 :   if (p <= (ulong)n) return 0;
    3399       56679 :   v = Fl_powu(2, factorial_lval(n, 2), p);
    3400       56757 :   for (k = 1;; k++)
    3401      142302 :   {
    3402      199059 :     long m = n >> (k-1), l, i;
    3403      199059 :     ulong a = 1;
    3404      199059 :     if (m <= 2) break;
    3405      142327 :     l = (1 + (n >> k)) | 1;
    3406             :     /* product of odd numbers in ]n / 2^k, 2 / 2^(k-1)] */
    3407      778458 :     for (i=l; i<=m; i+=2)
    3408      636144 :       a = Fl_mul(a, i, p);
    3409      142314 :     v = Fl_mul(v, k == 1? a: Fl_powu(a, k, p), p);
    3410             :   }
    3411       56732 :   return v;
    3412             : }
    3413             : 
    3414             : GEN
    3415         158 : factorial_Fp(long n, GEN p)
    3416             : {
    3417         158 :   pari_sp av = avma;
    3418             :   long k;
    3419         158 :   GEN v = Fp_powu(gen_2, factorial_lval(n, 2), p);
    3420         158 :   for (k = 1;; k++)
    3421         344 :   {
    3422         502 :     long m = n >> (k-1), l, i;
    3423         502 :     GEN a = gen_1;
    3424         502 :     if (m <= 2) break;
    3425         344 :     l = (1 + (n >> k)) | 1;
    3426             :     /* product of odd numbers in ]n / 2^k, 2 / 2^(k-1)] */
    3427         850 :     for (i=l; i<=m; i+=2)
    3428         506 :       a = Fp_mulu(a, i, p);
    3429         344 :     v = Fp_mul(v, k == 1? a: Fp_powu(a, k, p), p);
    3430         344 :     v = gerepileuptoint(av, v);
    3431             :   }
    3432         158 :   return v;
    3433             : }
    3434             : 
    3435             : /*******************************************************************/
    3436             : /**                      LUCAS & FIBONACCI                        **/
    3437             : /*******************************************************************/
    3438             : static void
    3439          56 : lucas(ulong n, GEN *a, GEN *b)
    3440             : {
    3441             :   GEN z, t, zt;
    3442          56 :   if (!n) { *a = gen_2; *b = gen_1; return; }
    3443          49 :   lucas(n >> 1, &z, &t); zt = mulii(z, t);
    3444          49 :   switch(n & 3) {
    3445          14 :     case  0: *a = subiu(sqri(z),2); *b = subiu(zt,1); break;
    3446          14 :     case  1: *a = subiu(zt,1);      *b = addiu(sqri(t),2); break;
    3447           7 :     case  2: *a = addiu(sqri(z),2); *b = addiu(zt,1); break;
    3448          14 :     case  3: *a = addiu(zt,1);      *b = subiu(sqri(t),2);
    3449             :   }
    3450             : }
    3451             : 
    3452             : GEN
    3453           7 : fibo(long n)
    3454             : {
    3455           7 :   pari_sp av = avma;
    3456             :   GEN a, b;
    3457           7 :   if (!n) return gen_0;
    3458           7 :   lucas((ulong)(labs(n)-1), &a, &b);
    3459           7 :   a = diviuexact(addii(shifti(a,1),b), 5);
    3460           7 :   if (n < 0 && !odd(n)) setsigne(a, -1);
    3461           7 :   return gerepileuptoint(av, a);
    3462             : }
    3463             : 
    3464             : /*******************************************************************/
    3465             : /*                      CONTINUED FRACTIONS                        */
    3466             : /*******************************************************************/
    3467             : static GEN
    3468     3136994 : icopy_lg(GEN x, long l)
    3469             : {
    3470     3136994 :   long lx = lgefint(x);
    3471             :   GEN y;
    3472             : 
    3473     3136994 :   if (lx >= l) return icopy(x);
    3474          49 :   y = cgeti(l); affii(x, y); return y;
    3475             : }
    3476             : 
    3477             : /* continued fraction of a/b. If y != NULL, stop when partial quotients
    3478             :  * differ from y */
    3479             : static GEN
    3480     3137344 : Qsfcont(GEN a, GEN b, GEN y, ulong k)
    3481             : {
    3482             :   GEN  z, c;
    3483     3137344 :   ulong i, l, ly = lgefint(b);
    3484             : 
    3485             :   /* times 1 / log2( (1+sqrt(5)) / 2 )  */
    3486     3137344 :   l = (ulong)(3 + bit_accuracy_mul(ly, 1.44042009041256));
    3487     3137344 :   if (k > 0 && k+1 > 0 && l > k+1) l = k+1; /* beware overflow */
    3488     3137344 :   if (l > LGBITS) l = LGBITS;
    3489             : 
    3490     3137344 :   z = cgetg(l,t_VEC);
    3491     3137344 :   l--;
    3492     3137344 :   if (y) {
    3493         350 :     pari_sp av = avma;
    3494         350 :     if (l >= (ulong)lg(y)) l = lg(y)-1;
    3495       25209 :     for (i = 1; i <= l; i++)
    3496             :     {
    3497       24985 :       GEN q = gel(y,i);
    3498       24985 :       gel(z,i) = q;
    3499       24985 :       c = b; if (!gequal1(q)) c = mulii(q, b);
    3500       24985 :       c = subii(a, c);
    3501       24985 :       if (signe(c) < 0)
    3502             :       { /* partial quotient too large */
    3503          96 :         c = addii(c, b);
    3504          96 :         if (signe(c) >= 0) i++; /* by 1 */
    3505          96 :         break;
    3506             :       }
    3507       24889 :       if (cmpii(c, b) >= 0)
    3508             :       { /* partial quotient too small */
    3509          30 :         c = subii(c, b);
    3510          30 :         if (cmpii(c, b) < 0) {
    3511             :           /* by 1. If next quotient is 1 in y, add 1 */
    3512          12 :           if (i < l && equali1(gel(y,i+1))) gel(z,i) = addiu(q,1);
    3513          12 :           i++;
    3514             :         }
    3515          30 :         break;
    3516             :       }
    3517       24859 :       if ((i & 0xff) == 0) gerepileall(av, 2, &b, &c);
    3518       24859 :       a = b; b = c;
    3519             :     }
    3520             :   } else {
    3521     3136994 :     a = icopy_lg(a, ly);
    3522     3136994 :     b = icopy(b);
    3523    24524282 :     for (i = 1; i <= l; i++)
    3524             :     {
    3525    24523964 :       gel(z,i) = truedvmdii(a,b,&c);
    3526    24523964 :       if (c == gen_0) { i++; break; }
    3527    21387288 :       affii(c, a); cgiv(c); c = a;
    3528    21387288 :       a = b; b = c;
    3529             :     }
    3530             :   }
    3531     3137344 :   i--;
    3532     3137344 :   if (i > 1 && gequal1(gel(z,i)))
    3533             :   {
    3534         101 :     cgiv(gel(z,i)); --i;
    3535         101 :     gel(z,i) = addui(1, gel(z,i)); /* unclean: leave old z[i] on stack */
    3536             :   }
    3537     3137344 :   setlg(z,i+1); return z;
    3538             : }
    3539             : 
    3540             : static GEN
    3541           0 : sersfcont(GEN a, GEN b, long k)
    3542             : {
    3543           0 :   long i, l = typ(a) == t_POL? lg(a): 3;
    3544             :   GEN y, c;
    3545           0 :   if (lg(b) > l) l = lg(b);
    3546           0 :   if (k > 0 && l > k+1) l = k+1;
    3547           0 :   y = cgetg(l,t_VEC);
    3548           0 :   for (i=1; i<l; i++)
    3549             :   {
    3550           0 :     gel(y,i) = poldivrem(a,b,&c);
    3551           0 :     if (gequal0(c)) { i++; break; }
    3552           0 :     a = b; b = c;
    3553             :   }
    3554           0 :   setlg(y, i); return y;
    3555             : }
    3556             : 
    3557             : GEN
    3558     3142293 : gboundcf(GEN x, long k)
    3559             : {
    3560             :   pari_sp av;
    3561     3142293 :   long tx = typ(x), e;
    3562             :   GEN y, a, b, c;
    3563             : 
    3564     3142293 :   if (k < 0) pari_err_DOMAIN("gboundcf","nmax","<",gen_0,stoi(k));
    3565     3142286 :   if (is_scalar_t(tx))
    3566             :   {
    3567     3142286 :     if (gequal0(x)) return mkvec(gen_0);
    3568     3142181 :     switch(tx)
    3569             :     {
    3570        5180 :       case t_INT: return mkveccopy(x);
    3571         357 :       case t_REAL:
    3572         357 :         av = avma;
    3573         357 :         c = mantissa_real(x,&e);
    3574         357 :         if (e < 0) pari_err_PREC("gboundcf");
    3575         350 :         y = int2n(e);
    3576         350 :         a = Qsfcont(c,y, NULL, k);
    3577         350 :         b = addsi(signe(x), c);
    3578         350 :         return gerepilecopy(av, Qsfcont(b,y, a, k));
    3579             : 
    3580     3136644 :       case t_FRAC:
    3581     3136644 :         av = avma;
    3582     3136644 :         return gerepileupto(av, Qsfcont(gel(x,1),gel(x,2), NULL, k));
    3583             :     }
    3584           0 :     pari_err_TYPE("gboundcf",x);
    3585             :   }
    3586             : 
    3587           0 :   switch(tx)
    3588             :   {
    3589           0 :     case t_POL: return mkveccopy(x);
    3590           0 :     case t_SER:
    3591           0 :       av = avma;
    3592           0 :       return gerepileupto(av, gboundcf(ser2rfrac_i(x), k));
    3593           0 :     case t_RFRAC:
    3594           0 :       av = avma;
    3595           0 :       return gerepilecopy(av, sersfcont(gel(x,1), gel(x,2), k));
    3596             :   }
    3597           0 :   pari_err_TYPE("gboundcf",x);
    3598             :   return NULL; /* LCOV_EXCL_LINE */
    3599             : }
    3600             : 
    3601             : static GEN
    3602          14 : sfcont2(GEN b, GEN x, long k)
    3603             : {
    3604          14 :   pari_sp av = avma;
    3605          14 :   long lb = lg(b), tx = typ(x), i;
    3606             :   GEN y,p1;
    3607             : 
    3608          14 :   if (k)
    3609             :   {
    3610           7 :     if (k >= lb) pari_err_DIM("contfrac [too few denominators]");
    3611           0 :     lb = k+1;
    3612             :   }
    3613           7 :   y = cgetg(lb,t_VEC);
    3614           7 :   if (lb==1) return y;
    3615           7 :   if (is_scalar_t(tx))
    3616             :   {
    3617           7 :     if (!is_intreal_t(tx) && tx != t_FRAC) pari_err_TYPE("sfcont2",x);
    3618             :   }
    3619           0 :   else if (tx == t_SER) x = ser2rfrac_i(x);
    3620             : 
    3621           7 :   if (!gequal1(gel(b,1))) x = gmul(gel(b,1),x);
    3622           7 :   for (i = 1;;)
    3623             :   {
    3624          35 :     if (tx == t_REAL)
    3625             :     {
    3626          35 :       long e = expo(x);
    3627          35 :       if (e > 0 && nbits2prec(e+1) > realprec(x)) break;
    3628          35 :       gel(y,i) = floorr(x);
    3629          35 :       p1 = subri(x, gel(y,i));
    3630             :     }
    3631             :     else
    3632             :     {
    3633           0 :       gel(y,i) = gfloor(x);
    3634           0 :       p1 = gsub(x, gel(y,i));
    3635             :     }
    3636          35 :     if (++i >= lb) break;
    3637          28 :     if (gequal0(p1)) break;
    3638          28 :     x = gdiv(gel(b,i),p1);
    3639             :   }
    3640           7 :   setlg(y,i);
    3641           7 :   return gerepilecopy(av,y);
    3642             : }
    3643             : 
    3644             : GEN
    3645         112 : gcf(GEN x) { return gboundcf(x,0); }
    3646             : GEN
    3647           0 : gcf2(GEN b, GEN x) { return contfrac0(x,b,0); }
    3648             : GEN
    3649          49 : contfrac0(GEN x, GEN b, long nmax)
    3650             : {
    3651             :   long tb;
    3652             : 
    3653          49 :   if (!b) return gboundcf(x,nmax);
    3654          28 :   tb = typ(b);
    3655          28 :   if (tb == t_INT) return gboundcf(x,itos(b));
    3656          21 :   if (! is_vec_t(tb)) pari_err_TYPE("contfrac0",b);
    3657          21 :   if (nmax < 0) pari_err_DOMAIN("contfrac","nmax","<",gen_0,stoi(nmax));
    3658          14 :   return sfcont2(b,x,nmax);
    3659             : }
    3660             : 
    3661             : GEN
    3662         252 : contfracpnqn(GEN x, long n)
    3663             : {
    3664         252 :   pari_sp av = avma;
    3665         252 :   long i, lx = lg(x);
    3666             :   GEN M,A,B, p0,p1, q0,q1;
    3667             : 
    3668         252 :   if (lx == 1)
    3669             :   {
    3670          28 :     if (! is_matvec_t(typ(x))) pari_err_TYPE("pnqn",x);
    3671          21 :     if (n >= 0) return cgetg(1,t_MAT);
    3672           7 :     return matid(2);
    3673             :   }
    3674         224 :   switch(typ(x))
    3675             :   {
    3676         182 :     case t_VEC: case t_COL: A = x; B = NULL; break;
    3677          42 :     case t_MAT:
    3678          42 :       switch(lgcols(x))
    3679             :       {
    3680           0 :         case 2: A = row(x,1); B = NULL; break;
    3681          35 :         case 3: A = row(x,2); B = row(x,1); break;
    3682           7 :         default: pari_err_DIM("pnqn [ nbrows != 1,2 ]");
    3683             :                  return NULL; /*LCOV_EXCL_LINE*/
    3684             :       }
    3685          35 :       break;
    3686           0 :     default: pari_err_TYPE("pnqn",x);
    3687             :       return NULL; /*LCOV_EXCL_LINE*/
    3688             :   }
    3689         217 :   p1 = gel(A,1);
    3690         217 :   q1 = B? gel(B,1): gen_1; /* p[0], q[0] */
    3691         217 :   if (n >= 0)
    3692             :   {
    3693         182 :     lx = minss(lx, n+2);
    3694         182 :     if (lx == 2) return gerepilecopy(av, mkmat(mkcol2(p1,q1)));
    3695             :   }
    3696          35 :   else if (lx == 2)
    3697           7 :     return gerepilecopy(av, mkmat2(mkcol2(p1,q1), mkcol2(gen_1,gen_0)));
    3698             :   /* lx >= 3 */
    3699         119 :   p0 = gen_1;
    3700         119 :   q0 = gen_0; /* p[-1], q[-1] */
    3701         119 :   M = cgetg(lx, t_MAT);
    3702         119 :   gel(M,1) = mkcol2(p1,q1);
    3703         399 :   for (i=2; i<lx; i++)
    3704             :   {
    3705         280 :     GEN a = gel(A,i), p2,q2;
    3706         280 :     if (B) {
    3707          84 :       GEN b = gel(B,i);
    3708          84 :       p0 = gmul(b,p0);
    3709          84 :       q0 = gmul(b,q0);
    3710             :     }
    3711         280 :     p2 = gadd(gmul(a,p1),p0); p0=p1; p1=p2;
    3712         280 :     q2 = gadd(gmul(a,q1),q0); q0=q1; q1=q2;
    3713         280 :     gel(M,i) = mkcol2(p1,q1);
    3714             :   }
    3715         119 :   if (n < 0) M = mkmat2(gel(M,lx-1), gel(M,lx-2));
    3716         119 :   return gerepilecopy(av, M);
    3717             : }
    3718             : GEN
    3719           0 : pnqn(GEN x) { return contfracpnqn(x,-1); }
    3720             : /* x = [a0, ..., an] from gboundcf, n >= 0;
    3721             :  * return [[p0, ..., pn], [q0,...,qn]] */
    3722             : GEN
    3723      894782 : ZV_allpnqn(GEN x)
    3724             : {
    3725      894782 :   long i, lx = lg(x);
    3726      894782 :   GEN p0, p1, q0, q1, p2, q2, P,Q, v = cgetg(3,t_VEC);
    3727             : 
    3728      894782 :   gel(v,1) = P = cgetg(lx, t_VEC);
    3729      894782 :   gel(v,2) = Q = cgetg(lx, t_VEC);
    3730      894782 :   p0 = gen_1; q0 = gen_0;
    3731      894782 :   gel(P, 1) = p1 = gel(x,1); gel(Q, 1) = q1 = gen_1;
    3732     3106138 :   for (i=2; i<lx; i++)
    3733             :   {
    3734     2211356 :     GEN a = gel(x,i);
    3735     2211356 :     gel(P, i) = p2 = addmulii(p0, a, p1); p0 = p1; p1 = p2;
    3736     2211356 :     gel(Q, i) = q2 = addmulii(q0, a, q1); q0 = q1; q1 = q2;
    3737             :   }
    3738      894782 :   return v;
    3739             : }
    3740             : 
    3741             : /* write Mod(x,N) as a/b, gcd(a,b) = 1, b <= B (no condition if B = NULL) */
    3742             : static GEN
    3743          42 : mod_to_frac(GEN x, GEN N, GEN B)
    3744             : {
    3745             :   GEN a, b, A;
    3746          42 :   if (B) A = divii(shifti(N, -1), B);
    3747             :   else
    3748             :   {
    3749          14 :     A = sqrti(shifti(N, -1));
    3750          14 :     B = A;
    3751             :   }
    3752          42 :   if (!Fp_ratlift(x, N, A,B,&a,&b) || !equali1( gcdii(a,b) )) return NULL;
    3753          28 :   return equali1(b)? a: mkfrac(a,b);
    3754             : }
    3755             : 
    3756             : static GEN
    3757         112 : mod_to_rfrac(GEN x, GEN N, long B)
    3758             : {
    3759             :   GEN a, b;
    3760         112 :   long A, d = degpol(N);
    3761         112 :   if (B >= 0) A = d-1 - B;
    3762             :   else
    3763             :   {
    3764          42 :     B = d >> 1;
    3765          42 :     A = odd(d)? B : B-1;
    3766             :   }
    3767         112 :   if (varn(N) != varn(x)) x = scalarpol(x, varn(N));
    3768         112 :   if (!RgXQ_ratlift(x, N, A, B, &a,&b) || degpol(RgX_gcd(a,b)) > 0) return NULL;
    3769          91 :   return gdiv(a,b);
    3770             : }
    3771             : 
    3772             : /* k > 0 t_INT, x a t_FRAC, returns the convergent a/b
    3773             :  * of the continued fraction of x with b <= k maximal */
    3774             : static GEN
    3775           7 : bestappr_frac(GEN x, GEN k)
    3776             : {
    3777             :   pari_sp av;
    3778             :   GEN p0, p1, p, q0, q1, q, a, y;
    3779             : 
    3780           7 :   if (cmpii(gel(x,2),k) <= 0) return gcopy(x);
    3781           0 :   av = avma; y = x;
    3782           0 :   p1 = gen_1; p0 = truedvmdii(gel(x,1), gel(x,2), &a); /* = floor(x) */
    3783           0 :   q1 = gen_0; q0 = gen_1;
    3784           0 :   x = mkfrac(a, gel(x,2)); /* = frac(x); now 0<= x < 1 */
    3785             :   for(;;)
    3786             :   {
    3787           0 :     x = ginv(x); /* > 1 */
    3788           0 :     a = typ(x)==t_INT? x: divii(gel(x,1), gel(x,2));
    3789           0 :     if (cmpii(a,k) > 0)
    3790             :     { /* next partial quotient will overflow limits */
    3791             :       GEN n, d;
    3792           0 :       a = divii(subii(k, q1), q0);
    3793           0 :       p = addii(mulii(a,p0), p1); p1=p0; p0=p;
    3794           0 :       q = addii(mulii(a,q0), q1); q1=q0; q0=q;
    3795             :       /* compare |y-p0/q0|, |y-p1/q1| */
    3796           0 :       n = gel(y,1);
    3797           0 :       d = gel(y,2);
    3798           0 :       if (abscmpii(mulii(q1, subii(mulii(q0,n), mulii(d,p0))),
    3799             :                    mulii(q0, subii(mulii(q1,n), mulii(d,p1)))) < 0)
    3800           0 :                    { p1 = p0; q1 = q0; }
    3801           0 :       break;
    3802             :     }
    3803           0 :     p = addii(mulii(a,p0), p1); p1=p0; p0=p;
    3804           0 :     q = addii(mulii(a,q0), q1); q1=q0; q0=q;
    3805             : 
    3806           0 :     if (cmpii(q0,k) > 0) break;
    3807           0 :     x = gsub(x,a); /* 0 <= x < 1 */
    3808           0 :     if (typ(x) == t_INT) { p1 = p0; q1 = q0; break; } /* x = 0 */
    3809             : 
    3810             :   }
    3811           0 :   return gerepileupto(av, gdiv(p1,q1));
    3812             : }
    3813             : /* k > 0 t_INT, x != 0 a t_REAL, returns the convergent a/b
    3814             :  * of the continued fraction of x with b <= k maximal */
    3815             : static GEN
    3816     1196052 : bestappr_real(GEN x, GEN k)
    3817             : {
    3818     1196052 :   pari_sp av = avma;
    3819     1196052 :   GEN kr, p0, p1, p, q0, q1, q, a, y = x;
    3820             : 
    3821     1196052 :   p1 = gen_1; a = p0 = floorr(x);
    3822     1195970 :   q1 = gen_0; q0 = gen_1;
    3823     1195970 :   x = subri(x,a); /* 0 <= x < 1 */
    3824     1196008 :   if (!signe(x)) { cgiv(x); return a; }
    3825     1082478 :   kr = itor(k, realprec(x));
    3826             :   for(;;)
    3827     1149692 :   {
    3828             :     long d;
    3829     2232192 :     x = invr(x); /* > 1 */
    3830     2232018 :     if (cmprr(x,kr) > 0)
    3831             :     { /* next partial quotient will overflow limits */
    3832     1060286 :       a = divii(subii(k, q1), q0);
    3833     1060318 :       p = addii(mulii(a,p0), p1); p1=p0; p0=p;
    3834     1060302 :       q = addii(mulii(a,q0), q1); q1=q0; q0=q;
    3835             :       /* compare |y-p0/q0|, |y-p1/q1| */
    3836     1060290 :       if (abscmprr(mulir(q1, subri(mulir(q0,y), p0)),
    3837             :                    mulir(q0, subri(mulir(q1,y), p1))) < 0)
    3838      127753 :                    { p1 = p0; q1 = q0; }
    3839     1060320 :       break;
    3840             :     }
    3841     1171829 :     d = nbits2prec(expo(x) + 1);
    3842     1171829 :     if (d > realprec(x)) { p1 = p0; q1 = q0; break; } /* original x was ~ 0 */
    3843             : 
    3844     1171639 :     a = truncr(x); /* truncr(x) will NOT raise e_PREC */
    3845     1171608 :     p = addii(mulii(a,p0), p1); p1=p0; p0=p;
    3846     1171630 :     q = addii(mulii(a,q0), q1); q1=q0; q0=q;
    3847             : 
    3848     1171630 :     if (cmpii(q0,k) > 0) break;
    3849     1165375 :     x = subri(x,a); /* 0 <= x < 1 */
    3850     1165382 :     if (!signe(x)) { p1 = p0; q1 = q0; break; }
    3851             :   }
    3852     1082455 :   if (signe(q1) < 0) { togglesign_safe(&p1); togglesign_safe(&q1); }
    3853     1082455 :   return gerepilecopy(av, equali1(q1)? p1: mkfrac(p1,q1));
    3854             : }
    3855             : 
    3856             : /* k t_INT or NULL */
    3857             : static GEN
    3858     2093791 : bestappr_Q(GEN x, GEN k)
    3859             : {
    3860     2093791 :   long lx, tx = typ(x), i;
    3861             :   GEN a, y;
    3862             : 
    3863     2093791 :   switch(tx)
    3864             :   {
    3865         154 :     case t_INT: return icopy(x);
    3866           7 :     case t_FRAC: return k? bestappr_frac(x, k): gcopy(x);
    3867     1427674 :     case t_REAL:
    3868     1427674 :       if (!signe(x)) return gen_0;
    3869             :       /* i <= e iff nbits2lg(e+1) > lg(x) iff floorr(x) fails */
    3870     1196036 :       i = bit_prec(x); if (i <= expo(x)) return NULL;
    3871     1196050 :       return bestappr_real(x, k? k: int2n(i));
    3872             : 
    3873          28 :     case t_INTMOD: {
    3874          28 :       pari_sp av = avma;
    3875          28 :       a = mod_to_frac(gel(x,2), gel(x,1), k); if (!a) return NULL;
    3876          21 :       return gerepilecopy(av, a);
    3877             :     }
    3878          14 :     case t_PADIC: {
    3879          14 :       pari_sp av = avma;
    3880          14 :       long v = valp(x);
    3881          14 :       a = mod_to_frac(gel(x,4), gel(x,3), k); if (!a) return NULL;
    3882           7 :       if (v) a = gmul(a, powis(gel(x,2), v));
    3883           7 :       return gerepilecopy(av, a);
    3884             :     }
    3885             : 
    3886        5453 :     case t_COMPLEX: {
    3887        5453 :       pari_sp av = avma;
    3888        5453 :       y = cgetg(3, t_COMPLEX);
    3889        5453 :       gel(y,2) = bestappr(gel(x,2), k);
    3890        5453 :       gel(y,1) = bestappr(gel(x,1), k);
    3891        5453 :       if (gequal0(gel(y,2))) return gerepileupto(av, gel(y,1));
    3892          91 :       return y;
    3893             :     }
    3894           0 :     case t_SER:
    3895           0 :       if (ser_isexactzero(x)) return gcopy(x);
    3896             :       /* fall through */
    3897             :     case t_POLMOD: case t_POL: case t_RFRAC:
    3898             :     case t_VEC: case t_COL: case t_MAT:
    3899      660461 :       y = cgetg_copy(x, &lx);
    3900      660507 :       if (lontyp[tx] == 1) i = 1; else { y[1] = x[1]; i = 2; }
    3901     2659002 :       for (; i<lx; i++)
    3902             :       {
    3903     1998517 :         a = bestappr_Q(gel(x,i),k); if (!a) return NULL;
    3904     1998495 :         gel(y,i) = a;
    3905             :       }
    3906      660485 :       if (tx == t_POL) return normalizepol(y);
    3907      660471 :       if (tx == t_SER) return normalizeser(y);
    3908      660471 :       return y;
    3909             :   }
    3910           0 :   pari_err_TYPE("bestappr_Q",x);
    3911             :   return NULL; /* LCOV_EXCL_LINE */
    3912             : }
    3913             : 
    3914             : static GEN
    3915          98 : bestappr_ser(GEN x, long B)
    3916             : {
    3917          98 :   long dN, v = valser(x), lx = lg(x);
    3918             :   GEN t;
    3919          98 :   x = normalizepol(ser2pol_i(x, lx));
    3920          98 :   dN = lx-2;
    3921          98 :   if (v > 0)
    3922             :   {
    3923          21 :     x = RgX_shift_shallow(x, v);
    3924          21 :     dN += v;
    3925             :   }
    3926          77 :   else if (v < 0)
    3927             :   {
    3928          14 :     if (B >= 0) B = maxss(B+v, 0);
    3929             :   }
    3930          98 :   t = mod_to_rfrac(x, pol_xn(dN, varn(x)), B);
    3931          98 :   if (!t) return NULL;
    3932          77 :   if (v < 0)
    3933             :   {
    3934             :     GEN a, b;
    3935             :     long vx;
    3936          14 :     if (typ(t) == t_POL) return RgX_mulXn(t, v);
    3937             :     /* t_RFRAC */
    3938          14 :     vx = varn(x);
    3939          14 :     a = gel(t,1);
    3940          14 :     b = gel(t,2);
    3941          14 :     v -= RgX_valrem(b, &b);
    3942          14 :     if (typ(a) == t_POL && varn(a) == vx) v += RgX_valrem(a, &a);
    3943          14 :     if (v < 0) b = RgX_shift_shallow(b, -v);
    3944           0 :     else if (v > 0) {
    3945           0 :       if (typ(a) != t_POL || varn(a) != vx) a = scalarpol_shallow(a, vx);
    3946           0 :       a = RgX_shift_shallow(a, v);
    3947             :     }
    3948          14 :     t = mkrfraccopy(a, b);
    3949             :   }
    3950          77 :   return t;
    3951             : }
    3952             : static GEN
    3953          42 : gc_empty(pari_sp av) { set_avma(av); return cgetg(1, t_VEC); }
    3954             : static GEN
    3955         112 : _gc_upto(pari_sp av, GEN x) { return x? gerepileupto(av, x): NULL; }
    3956             : 
    3957             : static GEN bestappr_RgX(GEN x, long B);
    3958             : /* B >= 0 or < 0 [omit condition on B].
    3959             :  * Look for coprime t_POL a,b, deg(b)<=B, such that a/b ~ x */
    3960             : static GEN
    3961         119 : bestappr_RgX(GEN x, long B)
    3962             : {
    3963             :   pari_sp av;
    3964         119 :   switch(typ(x))
    3965             :   {
    3966           0 :     case t_INT: case t_REAL: case t_INTMOD: case t_FRAC: case t_FFELT:
    3967             :     case t_COMPLEX: case t_PADIC: case t_QUAD: case t_POL:
    3968           0 :       return gcopy(x);
    3969          14 :     case t_RFRAC:
    3970          14 :       if (B < 0 || degpol(gel(x,2)) <= B) return gcopy(x);
    3971           7 :       av = avma; return _gc_upto(av, bestappr_ser(rfrac_to_ser_i(x, 2*B+1), B));
    3972          14 :     case t_POLMOD:
    3973          14 :       av = avma; return _gc_upto(av, mod_to_rfrac(gel(x,2), gel(x,1), B));
    3974          91 :     case t_SER:
    3975          91 :       av = avma; return _gc_upto(av, bestappr_ser(x, B));
    3976           0 :     case t_VEC: case t_COL: case t_MAT: {
    3977             :       long i, lx;
    3978           0 :       GEN y = cgetg_copy(x, &lx);
    3979           0 :       for (i = 1; i < lx; i++)
    3980             :       {
    3981           0 :         GEN t = bestappr_RgX(gel(x,i),B); if (!t) return NULL;
    3982           0 :         gel(y,i) = t;
    3983             :       }
    3984           0 :       return y;
    3985             :     }
    3986             :   }
    3987           0 :   pari_err_TYPE("bestappr_RgX",x);
    3988             :   return NULL; /* LCOV_EXCL_LINE */
    3989             : }
    3990             : 
    3991             : /* allow k = NULL: maximal accuracy */
    3992             : GEN
    3993       95270 : bestappr(GEN x, GEN k)
    3994             : {
    3995       95270 :   pari_sp av = avma;
    3996       95270 :   if (k) { /* replace by floor(k) */
    3997       94948 :     switch(typ(k))
    3998             :     {
    3999       33026 :       case t_INT:
    4000       33026 :         break;
    4001       61922 :       case t_REAL: case t_FRAC:
    4002       61922 :         k = floor_safe(k); /* left on stack for efficiency */
    4003       61925 :         if (!signe(k)) k = gen_1;
    4004       61925 :         break;
    4005           0 :       default:
    4006           0 :         pari_err_TYPE("bestappr [bound type]", k);
    4007           0 :         break;
    4008             :     }
    4009             :   }
    4010       95273 :   x = bestappr_Q(x, k);
    4011       95271 :   return x? x: gc_empty(av);
    4012             : }
    4013             : GEN
    4014         119 : bestapprPade(GEN x, long B)
    4015             : {
    4016         119 :   pari_sp av = avma;
    4017         119 :   GEN t = bestappr_RgX(x, B);
    4018         119 :   return t? t: gc_empty(av);
    4019             : }
    4020             : 
    4021             : static GEN
    4022          49 : serPade(GEN S, long p, long q)
    4023             : {
    4024          49 :   pari_sp av = avma;
    4025          49 :   long va, v, t = typ(S);
    4026          49 :   if (t!=t_SER && t!=t_POL && t!=t_RFRAC) pari_err_TYPE("bestapprPade", S);
    4027          49 :   va = gvar(S); v = gvaluation(S, pol_x(va));
    4028          49 :   if (p < 0) pari_err_DOMAIN("bestapprPade", "p", "<", gen_0, stoi(p));
    4029          49 :   if (q < 0) pari_err_DOMAIN("bestapprPade", "q", "<", gen_0, stoi(q));
    4030          49 :   if (v == LONG_MAX) return gc_empty(av);
    4031          42 :   S = gadd(S, zeroser(va, p + q + 1 + v));
    4032          42 :   return gerepileupto(av, bestapprPade(S, q));
    4033             : }
    4034             : 
    4035             : GEN
    4036         126 : bestapprPade0(GEN x, long p, long q)
    4037             : {
    4038          77 :   return (p >= 0 && q >= 0)? serPade(x, p, q)
    4039         203 :                            : bestapprPade(x, p >= 0? p: q);
    4040             : }

Generated by: LCOV version 1.16