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.18.1 lcov report (development 30005-fc14bb602a) Lines: 2084 2287 91.1 %
Date: 2025-02-18 09:22:46 Functions: 219 233 94.0 %
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      464140 : 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      464140 : 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      531469 : u_is_gener_expo(ulong p, GEN L0)
      54             : {
      55      531469 :   const ulong q = p >> 1;
      56             :   long i;
      57             :   GEN L;
      58      531469 :   if (!L0) L0 = u_odd_prime_divisors(q);
      59      531471 :   L = cgetg_copy(L0,&i);
      60     1150448 :   while (--i) L[i] = q / uel(L0,i);
      61      531471 :   return L;
      62             : }
      63             : 
      64             : int
      65     1628850 : is_gener_Fl(ulong x, ulong p, ulong p_1, GEN L)
      66             : {
      67             :   long i;
      68     1628850 :   if (krouu(x, p) >= 0) return 0;
      69     1375560 :   for (i=lg(L)-1; i; i--)
      70             :   {
      71      839352 :     ulong t = Fl_powu(x, uel(L,i), p);
      72      839352 :     if (t == p_1 || t == 1) return 0;
      73             :   }
      74      536208 :   return 1;
      75             : }
      76             : /* assume p prime */
      77             : ulong
      78     1051268 : pgener_Fl_local(ulong p, GEN L0)
      79             : {
      80     1051268 :   const pari_sp av = avma;
      81     1051268 :   const ulong p_1 = p-1;
      82             :   long x;
      83             :   GEN L;
      84     1051268 :   if (p <= 19) switch(p)
      85             :   { /* quick trivial cases */
      86          63 :     case 2:  return 1;
      87      116195 :     case 7:
      88      116195 :     case 17: return 3;
      89      403586 :     default: return 2;
      90             :   }
      91      531424 :   L = u_is_gener_expo(p,L0);
      92     1621090 :   for (x = 2;; x++)
      93     1621090 :     if (is_gener_Fl(x,p,p_1,L)) return gc_ulong(av, x);
      94             : }
      95             : ulong
      96      575230 : 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       13742 : is_gener_Fp(GEN x, GEN p, GEN p_1, GEN L)
     102             : {
     103       13742 :   long i, t = lgefint(x)==3? kroui(x[2], p): kronecker(x, p);
     104       13742 :   if (t >= 0) return 0;
     105       21739 :   for (i = lg(L)-1; i; i--)
     106             :   {
     107       14237 :     GEN t = Fp_pow(x, gel(L,i), p);
     108       14237 :     if (equalii(t, p_1) || equali1(t)) return 0;
     109             :   }
     110        7502 :   return 1;
     111             : }
     112             : 
     113             : /* assume p prime, return a generator of all L[i]-Sylows in F_p^*. */
     114             : GEN
     115      358163 : pgener_Fp_local(GEN p, GEN L0)
     116             : {
     117      358163 :   pari_sp av0 = avma;
     118             :   GEN x, p_1, L;
     119      358163 :   if (lgefint(p) == 3)
     120             :   {
     121             :     ulong z;
     122      353260 :     if (p[2] == 2) return gen_1;
     123      258273 :     if (L0) L0 = ZV_to_nv(L0);
     124      258272 :     z = pgener_Fl_local(uel(p,2), L0);
     125      258306 :     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       44247 : pgener_Fp(GEN p) { return pgener_Fp_local(p, NULL); }
     135             : 
     136             : ulong
     137      205710 : pgener_Zl(ulong p)
     138             : {
     139      205710 :   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      205710 :   if (p == 40487) return 10;
     142             : #ifndef LONG_IS_64BIT
     143       29808 :   return pgener_Fl(p);
     144             : #else
     145      175902 :   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      170954 : pgener_Zp(GEN p)
     162             : {
     163      170954 :   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      215568 : rootsof1_Fl(ulong n, ulong p)
     251             : {
     252      215568 :   pari_sp av = avma;
     253      215568 :   GEN L = u_odd_prime_divisors(n); /* 2 implicit in pgener_Fl_local */
     254      215568 :   ulong z = pgener_Fl_local(p, L);
     255      215568 :   z = Fl_powu(z, (p-1) / n, p); /* prim. n-th root of 1 */
     256      215568 :   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   321343056 : ome(long t)
     364             : {
     365   321343056 :   switch(t & 7)
     366             :   {
     367   182293421 :     case 3:
     368   182293421 :     case 5: return 1;
     369   139049635 :     default: return 0;
     370             :   }
     371             : }
     372             : /* t a t_INT, is t = 3,5 mod 8 ? */
     373             : static int
     374     5605703 : gome(GEN t)
     375     5605703 : { return signe(t)? ome( mod2BIL(t) ): 0; }
     376             : 
     377             : /* assume y odd, return kronecker(x,y) * s */
     378             : static long
     379   228161078 : krouu_s(ulong x, ulong y, long s)
     380             : {
     381   228161078 :   ulong x1 = x, y1 = y, z;
     382  1033028537 :   while (x1)
     383             :   {
     384   804864550 :     long r = vals(x1);
     385   804962566 :     if (r)
     386             :     {
     387   427948271 :       if (odd(r) && ome(y1)) s = -s;
     388   427853164 :       x1 >>= r;
     389             :     }
     390   804867459 :     if (x1 & y1 & 2) s = -s;
     391   804867459 :     z = y1 % x1; y1 = x1; x1 = z;
     392             :   }
     393   228163987 :   return (y1 == 1)? s: 0;
     394             : }
     395             : 
     396             : long
     397    11963112 : kronecker(GEN x, GEN y)
     398             : {
     399    11963112 :   pari_sp av = avma;
     400    11963112 :   long s = 1, r;
     401             :   ulong xu;
     402             : 
     403    11963112 :   if (typ(x) != t_INT) pari_err_TYPE("kronecker",x);
     404    11963112 :   if (typ(y) != t_INT) pari_err_TYPE("kronecker",y);
     405    11963112 :   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    11962979 :   r = vali(y);
     411    11962975 :   if (r)
     412             :   {
     413     1348912 :     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    10935774 :   x = modii(x,y);
     418    13351449 :   while (lgefint(x) > 3) /* x < y */
     419             :   {
     420             :     GEN z;
     421     2415710 :     r = vali(x);
     422     2415363 :     if (r)
     423             :     {
     424     1318641 :       if (odd(r) && gome(y)) s = -s;
     425     1318564 :       x = shifti(x,-r);
     426             :     }
     427             :     /* x=3 mod 4 && y=3 mod 4 ? (both are odd here) */
     428     2414674 :     if (mod2BIL(x) & mod2BIL(y) & 2) s = -s;
     429     2414061 :     z = remii(y,x); y = x; x = z;
     430     2415634 :     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    10935739 :   xu = itou(x);
     437    10935734 :   if (!xu) return is_pm1(y)? s: 0;
     438    10837649 :   r = vals(xu);
     439    10837647 :   if (r)
     440             :   {
     441     5754927 :     if (odd(r) && gome(y)) s = -s;
     442     5754927 :     xu >>= r;
     443             :   }
     444             :   /* x=3 mod 4 && y=3 mod 4 ? (both are odd here) */
     445    10837647 :   if (xu & mod2BIL(y) & 2) s = -s;
     446    10837649 :   return gc_long(av, krouu_s(umodiu(y,xu), xu, s));
     447             : }
     448             : 
     449             : long
     450       39753 : krois(GEN x, long y)
     451             : {
     452             :   ulong yu;
     453       39753 :   long s = 1;
     454             : 
     455       39753 :   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       39725 :     yu = (ulong)y;
     462       39725 :   if (!odd(yu))
     463             :   {
     464             :     long r;
     465       18417 :     if (!mpodd(x)) return 0;
     466       12467 :     r = vals(yu); yu >>= r;
     467       12467 :     if (odd(r) && gome(x)) s = -s;
     468             :   }
     469       33775 :   return krouu_s(umodiu(x, yu), yu, s);
     470             : }
     471             : /* assume y != 0 */
     472             : long
     473    27620405 : kroiu(GEN x, ulong y)
     474             : {
     475             :   long r;
     476    27620405 :   if (odd(y)) return krouu_s(umodiu(x,y), y, 1);
     477      302998 :   if (!mpodd(x)) return 0;
     478      208344 :   r = vals(y); y >>= r;
     479      208345 :   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      178046 : krouodd(ulong x, GEN y, long s)
     485             : {
     486             :   long r;
     487      178046 :   if (lgefint(y) == 3) return krouu_s(x, y[2], s);
     488       27958 :   if (!x) return 0; /* y != 1 */
     489       27958 :   r = vals(x);
     490       27958 :   if (r)
     491             :   {
     492       14499 :     if (odd(r) && gome(y)) s = -s;
     493       14499 :     x >>= r;
     494             :   }
     495             :   /* x=3 mod 4 && y=3 mod 4 ? (both are odd here) */
     496       27958 :   if (x & mod2BIL(y) & 2) s = -s;
     497       27958 :   return krouu_s(umodiu(y,x), x, s);
     498             : }
     499             : 
     500             : long
     501      143221 : krosi(long x, GEN y)
     502             : {
     503      143221 :   const pari_sp av = avma;
     504      143221 :   long s = 1, r;
     505      143221 :   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      143221 :   r = vali(y);
     511      143221 :   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      143221 :   if (x < 0) { x = -x; if (mod4(y) == 3) s = -s; }
     518      143221 :   return gc_long(av, krouodd((ulong)x, y, s));
     519             : }
     520             : 
     521             : long
     522       34825 : kroui(ulong x, GEN y)
     523             : {
     524       34825 :   const pari_sp av = avma;
     525       34825 :   long s = 1, r;
     526       34825 :   switch (signe(y))
     527             :   {
     528           0 :     case -1: y = negi(y); break;
     529           0 :     case 0: return x==1UL;
     530             :   }
     531       34825 :   r = vali(y);
     532       34825 :   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       34825 :   return gc_long(av, krouodd(x, y, s));
     539             : }
     540             : 
     541             : long
     542    97777671 : kross(long x, long y)
     543             : {
     544             :   ulong yu;
     545    97777671 :   long s = 1;
     546             : 
     547    97777671 :   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    97708728 :     yu = (ulong)y;
     554    97777643 :   if (!odd(yu))
     555             :   {
     556             :     long r;
     557    23582767 :     if (!odd(x)) return 0;
     558    16659929 :     r = vals(yu); yu >>= r;
     559    16659929 :     if (odd(r) && ome(x)) s = -s;
     560             :   }
     561    90854805 :   x %= (long)yu; if (x < 0) x += yu;
     562    90854805 :   return krouu_s((ulong)x, yu, s);
     563             : }
     564             : 
     565             : long
     566    98746808 : krouu(ulong x, ulong y)
     567             : {
     568             :   long r;
     569    98746808 :   if (odd(y)) return krouu_s(x, y, 1);
     570       16953 :   if (!odd(x)) return 0;
     571       17037 :   r = vals(y); y >>= r;
     572       17037 :   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      140826 : hilbertii(GEN x, GEN y, GEN p)
     589             : {
     590             :   pari_sp av;
     591             :   long oddvx, oddvy, z;
     592             : 
     593      140826 :   if (!p) return mphilbertoo(x,y);
     594      133518 :   if (is_pm1(p) || signe(p) < 0) pari_err_PRIME("hilbertii",p);
     595      133518 :   if (!signe(x) || !signe(y)) return 0;
     596      133497 :   av = avma;
     597      133497 :   oddvx = odd(Z_pvalrem(x,p,&x));
     598      133497 :   oddvy = odd(Z_pvalrem(y,p,&y));
     599             :   /* x, y are p-units, compute hilbert(x * p^oddvx, y * p^oddvy, p) */
     600      133497 :   if (absequaliu(p, 2))
     601             :   {
     602       12355 :     z = (Mod4(x) == 3 && Mod4(y) == 3)? -1: 1;
     603       12355 :     if (oddvx && gome(y)) z = -z;
     604       12355 :     if (oddvy && gome(x)) z = -z;
     605             :   }
     606             :   else
     607             :   {
     608      121142 :     z = (oddvx && oddvy && mod4(p) == 3)? -1: 1;
     609      121142 :     if (oddvx && kronecker(y,p) < 0) z = -z;
     610      121142 :     if (oddvy && kronecker(x,p) < 0) z = -z;
     611             :   }
     612      133497 :   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 = padic_p(x), u = padic_u(x);
     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(u)) err_prec();
     659          70 :   return odd(valp(x))? mulii(p,u): u;
     660             : }
     661             : 
     662             : long
     663       62314 : hilbert(GEN x, GEN y, GEN p)
     664             : {
     665       62314 :   pari_sp av = avma;
     666       62314 :   long tx = typ(x), ty = typ(y);
     667             : 
     668       62314 :   if (p && typ(p) != t_INT) pari_err_TYPE("hilbert",p);
     669       62314 :   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       62237 :   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       62223 :   if (tx == t_INTMOD) { x = lift_intmod(x, &p); tx = t_INT; }
     692       62020 :   if (ty == t_INTMOD) { y = lift_intmod(y, &p); ty = t_INT; }
     693             : 
     694       61964 :   if (tx == t_PADIC) { x = lift_padic(x, &p); tx = t_INT; }
     695       61901 :   if (ty == t_PADIC) { y = lift_padic(y, &p); ty = t_INT; }
     696             : 
     697       61824 :   if (tx == t_FRAC) { tx = t_INT; x = p? mulii(gel(x,1),gel(x,2)): gel(x,1); }
     698       61824 :   if (ty == t_FRAC) { ty = t_INT; y = p? mulii(gel(y,1),gel(y,2)): gel(y,1); }
     699             : 
     700       61824 :   if (tx != t_INT || ty != t_INT) pari_err_TYPE2("hilbert",x,y);
     701       61824 :   if (p && !signe(p)) p = NULL;
     702       61824 :   return gc_long(av, hilbertii(x,y,p));
     703             : }
     704             : 
     705             : /*******************************************************************/
     706             : /*                       SQUARE ROOT MODULO p                      */
     707             : /*******************************************************************/
     708             : static void
     709     2257030 : checkp(ulong q, ulong p)
     710     2257030 : { 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    11668502 : nonsquare1_Fl(ulong p)
     714             : {
     715             :   forprime_t S;
     716             :   ulong q;
     717    11668502 :   if ((p & 7UL) != 1) return 2UL;
     718     4360728 :   q = p % 3; if (q == 2) return 3UL;
     719     1403212 :   checkp(q, p);
     720     1410358 :   q = p % 5; if (q == 2 || q == 3) return 5UL;
     721      522816 :   checkp(q, p);
     722      522810 :   q = p % 7; if (q != 4 && q >= 3) return 7UL;
     723      195597 :   checkp(q, p);
     724             :   /* log^2(2^64) < 1968 is enough under GRH (and p^(1/4)log(p) without it)*/
     725      195648 :   u_forprime_init(&S, 11, 1967);
     726      324077 :   while ((q = u_forprime_next(&S)))
     727             :   {
     728      324074 :     if (krouu(q, p) < 0) return q;
     729      128426 :     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      177391 : Fl_2gener_pre(ulong p, ulong pi)
     742             : {
     743      177391 :   ulong p1 = p-1;
     744      177391 :   long e = vals(p1);
     745      177372 :   if (e == 1) return p1;
     746       65805 :   return Fl_powu_pre(nonsquare1_Fl(p), p1 >> e, p, pi);
     747             : }
     748             : 
     749             : ulong
     750       67676 : Fl_2gener_pre_i(ulong  ns, ulong p, ulong pi)
     751             : {
     752       67676 :   ulong p1 = p-1;
     753       67676 :   long e = vals(p1);
     754       67676 :   if (e == 1) return p1;
     755       25875 :   return Fl_powu_pre(ns, p1 >> e, p, pi);
     756             : }
     757             : 
     758             : static ulong
     759    12454837 : Fl_sqrt_i(ulong a, ulong y, ulong p)
     760             : {
     761             :   long i, e, k;
     762             :   ulong p1, q, v, w;
     763             : 
     764    12454837 :   if (!a) return 0;
     765    11175371 :   p1 = p - 1; e = vals(p1);
     766    11176071 :   if (e == 0) /* p = 2 */
     767             :   {
     768      650263 :     if (p != 2) pari_err_PRIME("Fl_sqrt [modulus]",utoi(p));
     769      651335 :     return ((a & 1) == 0)? 0: 1;
     770             :   }
     771    10525808 :   if (e == 1)
     772             :   {
     773     5007052 :     v = Fl_powu(a, (p+1) >> 2, p);
     774     5007136 :     if (Fl_sqr(v, p) != a) return ~0UL;
     775     5002274 :     p1 = p - v; if (v > p1) v = p1;
     776     5002274 :     return v;
     777             :   }
     778     5518756 :   q = p1 >> e; /* q = (p-1)/2^oo is odd */
     779     5518756 :   p1 = Fl_powu(a, q >> 1, p); /* a ^ [(q-1)/2] */
     780     5518797 :   if (!p1) return 0;
     781     5518797 :   v = Fl_mul(a, p1, p);
     782     5518803 :   w = Fl_mul(v, p1, p);
     783     5518823 :   if (!y) y = Fl_powu(nonsquare1_Fl(p), q, p);
     784     9371961 :   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     3855048 :     p1 = Fl_sqr(w, p);
     788     6385251 :     for (k=1; p1 != 1 && k < e; k++) p1 = Fl_sqr(p1, p);
     789     3855052 :     if (k == e) return ~0UL;
     790             :     /* w ^ (2^k) = 1 --> w = y ^ (u * 2^(e-k)), u odd */
     791     3852975 :     p1 = y;
     792     5130449 :     for (i=1; i < e-k; i++) p1 = Fl_sqr(p1, p);
     793     3852978 :     y = Fl_sqr(p1, p); e = k;
     794     3853001 :     w = Fl_mul(y, w, p);
     795     3853015 :     v = Fl_mul(v, p1, p);
     796             :   }
     797     5516913 :   p1 = p - v; if (v > p1) v = p1;
     798     5516913 :   return v;
     799             : }
     800             : 
     801             : /* Tonelli-Shanks. Assume p is prime and (a,p) != -1. Allow pi = 0 */
     802             : ulong
     803    33734718 : 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    33734718 :   if (!pi) return Fl_sqrt_i(a, y, p);
     809    21279892 :   if (!a) return 0;
     810    21153896 :   p1 = p - 1; e = vals(p1);
     811    21155280 :   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    21170858 :   if (e == 1)
     817             :   {
     818    15034903 :     v = Fl_powu_pre(a, (p+1) >> 2, p, pi);
     819    15025283 :     if (Fl_sqr_pre(v, p, pi) != a) return ~0UL;
     820    15031427 :     p1 = p - v; if (v > p1) v = p1;
     821    15031427 :     return v;
     822             :   }
     823     6135955 :   q = p1 >> e; /* q = (p-1)/2^oo is odd */
     824     6135955 :   p1 = Fl_powu_pre(a, q >> 1, p, pi); /* a ^ [(q-1)/2] */
     825     6134288 :   if (!p1) return 0;
     826     6134288 :   v = Fl_mul_pre(a, p1, p, pi);
     827     6134844 :   w = Fl_mul_pre(v, p1, p, pi);
     828     6133065 :   if (!y) y = Fl_powu_pre(nonsquare1_Fl(p), q, p, pi);
     829    11671044 :   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     5536240 :     p1 = Fl_sqr_pre(w,p,pi);
     833    10355499 :     for (k=1; p1 != 1 && k < e; k++) p1 = Fl_sqr_pre(p1,p,pi);
     834     5536725 :     if (k == e) return ~0UL;
     835             :     /* w ^ (2^k) = 1 --> w = y ^ (u * 2^(e-k)), u odd */
     836     5536633 :     p1 = y;
     837     7272568 :     for (i=1; i < e-k; i++) p1 = Fl_sqr_pre(p1, p, pi);
     838     5536661 :     y = Fl_sqr_pre(p1, p, pi); e = k;
     839     5538192 :     w = Fl_mul_pre(y, w, p, pi);
     840     5536483 :     v = Fl_mul_pre(v, p1, p, pi);
     841             :   }
     842     6134804 :   p1 = p - v; if (v > p1) v = p1;
     843     6134804 :   return v;
     844             : }
     845             : 
     846             : ulong
     847    12330298 : Fl_sqrt(ulong a, ulong p)
     848    12330298 : { ulong pi = (p & HIGHMASK)? get_Fl_red(p): 0; return Fl_sqrt_pre_i(a, 0, p, pi); }
     849             : 
     850             : ulong
     851    21212668 : Fl_sqrt_pre(ulong a, ulong p, ulong pi)
     852    21212668 : { return Fl_sqrt_pre_i(a, 0, p, pi); }
     853             : 
     854             : /* allow pi = 0 */
     855             : static ulong
     856      141882 : Fl_lgener_pre_all(ulong l, long e, ulong r, ulong p, ulong pi, ulong *pt_m)
     857             : {
     858      141882 :   ulong x, y, m, le1 = upowuu(l, e-1);
     859      141882 :   for (x = 2; ; x++)
     860             :   {
     861      173015 :     y = Fl_powu_pre(x, r, p, pi);
     862      173014 :     if (y==1) continue;
     863      155669 :     m = Fl_powu_pre(y, le1, p, pi);
     864      155671 :     if (m != 1) break;
     865             :   }
     866      141883 :   *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      226140 : 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      226140 :   if (a==0) return a;
     879      226134 :   u2 = Fl_inv(l%r, r);
     880      226134 :   v = Fl_powu_pre(a, u2, p, pi);
     881      226131 :   w = Fl_powu_pre(v, l, p, pi);
     882      226129 :   w = pi? Fl_mul_pre(w, Fl_inv(a, p), p, pi): Fl_div(w, a, p);
     883      226119 :   if (w==1) return v;
     884      139226 :   if (y==0) y = Fl_lgener_pre_all(l, e, r, p, pi, &m);
     885      164784 :   while (w!=1)
     886             :   {
     887      144666 :     ulong k = 0;
     888      144666 :     p1 = w;
     889             :     do
     890             :     {
     891      188297 :       z = p1; p1 = Fl_powu_pre(p1, l, p, pi);
     892      188298 :       if (++k == e) return ULONG_MAX;
     893       69188 :     } while (p1!=1);
     894       25557 :     dl = Fl_log_pre(z, m, l, p, pi);
     895       25557 :     dl = Fl_neg(dl, l);
     896       25557 :     p1 = Fl_powu_pre(y,dl*upowuu(l,e-k-1),p,pi);
     897       25557 :     m = Fl_powu_pre(m, dl, p, pi);
     898       25557 :     e = k;
     899       25557 :     v = pi? Fl_mul_pre(p1,v,p,pi): Fl_mul(p1,v,p);
     900       25557 :     y = Fl_powu_pre(p1,l,p,pi);
     901       25557 :     w = pi? Fl_mul_pre(y,w,p,pi): Fl_mul(y,w,p);
     902             :   }
     903       20118 :   return v;
     904             : }
     905             : 
     906             : /* allow pi = 0 */
     907             : static ulong
     908      223258 : Fl_sqrtl_i(ulong a, ulong l, ulong p, ulong pi, ulong y, ulong m)
     909             : {
     910      223258 :   ulong r, e = u_lvalrem(p-1, l, &r);
     911      223258 :   return Fl_sqrtl_raw(a, l, e, r, p, pi, y, m);
     912             : }
     913             : /* allow pi = 0 */
     914             : ulong
     915      223258 : Fl_sqrtl_pre(ulong a, ulong l, ulong p, ulong pi)
     916      223258 : { 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      232483 : Fl_sqrtn_pre(ulong a, long n, ulong p, ulong pi, ulong *zetan)
     926             : {
     927      232483 :   ulong m, q = p-1, z;
     928      232483 :   ulong nn = n >= 0 ? (ulong)n: -(ulong)n;
     929      232483 :   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      116094 :   if (n==1)
     936             :   {
     937         420 :     if (zetan) *zetan = 1;
     938         420 :     return n < 0? Fl_inv(a,p): a;
     939             :   }
     940      115674 :   if (n==2)
     941             :   {
     942       42082 :     if (zetan) *zetan = p-1;
     943       42082 :     return Fl_sqrt_pre_i(a, 0, p, pi);
     944             :   }
     945       73592 :   if (a == 1 && !zetan) return a;
     946       44136 :   m = ugcd(nn, q);
     947       44136 :   z = 1;
     948       44136 :   if (m!=1)
     949             :   {
     950        2668 :     GEN F = factoru(m);
     951             :     long i, j, e;
     952             :     ulong r, zeta, y, l;
     953        5733 :     for (i = nbrows(F); i; i--)
     954             :     {
     955        3114 :       l = ucoeff(F,i,1);
     956        3114 :       j = ucoeff(F,i,2);
     957        3114 :       e = u_lvalrem(q,l, &r);
     958        3114 :       y = Fl_lgener_pre_all(l, e, r, p, pi, &zeta);
     959        3114 :       if (zetan)
     960             :       {
     961        1354 :         ulong Y = Fl_powu_pre(y, upowuu(l,e-j), p, pi);
     962        1354 :         z = pi? Fl_mul_pre(z, Y, p, pi): Fl_mul(z, Y, p);
     963             :       }
     964        3114 :       if (a!=1)
     965             :         do
     966             :         {
     967        2880 :           a = Fl_sqrtl_raw(a, l, e, r, p, pi, y, zeta);
     968        2866 :           if (a==ULONG_MAX) return ULONG_MAX;
     969        2831 :         } while (--j);
     970             :     }
     971             :   }
     972       44087 :   if (m != nn)
     973             :   {
     974       41489 :     ulong qm = q/m, nm = (nn/m) % qm;
     975       41489 :     a = Fl_powu_pre(a, Fl_inv(nm, qm), p, pi);
     976             :   }
     977       44087 :   if (n < 0) a = Fl_inv(a, p);
     978       44087 :   if (zetan) *zetan = z;
     979       44087 :   return a;
     980             : }
     981             : 
     982             : ulong
     983      232483 : Fl_sqrtn(ulong a, long n, ulong p, ulong *zetan)
     984             : {
     985      232483 :   ulong pi = (p & HIGHMASK)? get_Fl_red(p): 0;
     986      232483 :   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        5914 : Fp_2gener_all(GEN q, GEN p)
    1056             : {
    1057             :   long k;
    1058        5914 :   for (k = 2;; k++)
    1059       11627 :   {
    1060       17541 :     long i = kroui(k, p);
    1061       17541 :     if (i < 0) return Fp_pow(utoipos(k), q, p);
    1062       11627 :     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       19821 : Fp_2gener_i(GEN ns, GEN p)
    1078             : {
    1079       19821 :   GEN q = subiu(p,1);
    1080       19821 :   long e = vali(q);
    1081       19821 :   if (e == 1) return q;
    1082       18567 :   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       19578 : Fp_sqrts(long a, GEN p)
    1133             : {
    1134       19578 :   long v = vals(a)>>1;
    1135       19578 :   GEN r = gen_0;
    1136       19578 :   a >>= v << 1;
    1137       19578 :   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        2214 :     case 5: case 13: case 17: case 21: case 29: case 33:
    1178             :     case -7: case -11: case -15: case -19: case -23:
    1179        2214 :       if (umodiu(p,labs(a))==1)
    1180         334 :         r = Fp_gausssum(a,p);
    1181             :       else
    1182        1880 :         return gen_0;
    1183         334 :       break;
    1184       15422 :     default:
    1185       15422 :       return gen_0;
    1186             :   }
    1187        2276 :   return remii(shifti(r, v), p);
    1188             : }
    1189             : 
    1190             : static GEN
    1191       78149 : Fp_sqrt_ii(GEN a, GEN y, GEN p)
    1192             : {
    1193       78149 :   pari_sp av = avma;
    1194       78149 :   GEN  q, v, w, p1 = subiu(p,1);
    1195       78151 :   long i, k, e = vali(p1), as;
    1196             : 
    1197             :   /* direct formulas more efficient */
    1198       78152 :   if (e == 0) pari_err_PRIME("Fp_sqrt [modulus]",p); /* p != 2 */
    1199       78152 :   if (e == 1)
    1200             :   {
    1201       19152 :     q = addiu(shifti(p1,-2),1); /* (p+1) / 4 */
    1202       19150 :     v = Fp_pow(a, q, p);
    1203             :     /* must check equality in case (a/p) = -1 or p not prime */
    1204       19158 :     av = avma; e = equalii(Fp_sqr(v,p), a); set_avma(av);
    1205       19156 :     return e? v: NULL;
    1206             :   }
    1207       59000 :   as = itos_or_0(a);
    1208       59002 :   if (!as) as = itos_or_0(subii(a,p));
    1209       59008 :   if (as)
    1210             :   {
    1211       19577 :     GEN res = Fp_sqrts(as, p);
    1212       19578 :     if (!res) return gc_NULL(av);
    1213       19578 :     if (signe(res)) return gerepileupto(av, res);
    1214             :   }
    1215       56733 :   if (e == 2)
    1216             :   { /* Atkin's formula */
    1217       17892 :     GEN I, a2 = shifti(a,1);
    1218       17892 :     if (cmpii(a2,p) >= 0) a2 = subii(a2,p);
    1219       17892 :     q = shifti(p1, -3); /* (p-5)/8 */
    1220       17892 :     v = Fp_pow(a2, q, p);
    1221       17893 :     I = Fp_mul(a2, Fp_sqr(v,p), p); /* I^2 = -1 */
    1222       17892 :     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       17893 :     av = avma; e = equalii(Fp_sqr(v,p), a); set_avma(av);
    1225       17892 :     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       38841 :   if (e*(e-1) > 20 + 8 * expi(p)) return sqrt_Cipolla(a,p);
    1230             :   /* Tonelli-Shanks */
    1231       38841 :   av = avma; q = shifti(p1,-e); /* q = (p-1)/2^oo is odd */
    1232       38840 :   if (!y)
    1233             :   {
    1234        2722 :     y = Fp_2gener_all(q, p);
    1235        2722 :     if (!y) pari_err_PRIME("Fp_sqrt [modulus]",p);
    1236             :   }
    1237       38840 :   p1 = Fp_pow(a, shifti(q,-1), p); /* a ^ (q-1)/2 */
    1238       38842 :   v = Fp_mul(a, p1, p);
    1239       38842 :   w = Fp_mul(v, p1, p);
    1240       92631 :   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       53832 :     p1 = Fp_sqr(w,p);
    1244      110448 :     for (k=1; !equali1(p1) && k < e; k++) p1 = Fp_sqr(p1,p);
    1245       53832 :     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       53789 :     p1 = y;
    1248       76543 :     for (i=1; i < e-k; i++) p1 = Fp_sqr(p1,p);
    1249       53789 :     y = Fp_sqr(p1, p); e = k;
    1250       53788 :     w = Fp_mul(y, w, p);
    1251       53789 :     v = Fp_mul(v, p1, p);
    1252       53789 :     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       38799 :   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     4891578 : Fp_sqrt_i(GEN a, GEN y, GEN p)
    1265             : {
    1266     4891578 :   pari_sp av = avma, av2;
    1267             :   GEN q;
    1268             : 
    1269     4891578 :   if (lgefint(p) == 3)
    1270             :   {
    1271     4813316 :     ulong pp = uel(p,2), u = umodiu(a, pp);
    1272     4813324 :     if (!u) return gen_0;
    1273     3600610 :     u = Fl_sqrt(u, pp);
    1274     3600729 :     return (u == ~0UL)? NULL: utoipos(u);
    1275             :   }
    1276       78262 :   a = modii(a, p); if (!signe(a)) return gen_0;
    1277       78149 :   a = Fp_sqrt_ii(a, y, p); if (!a) return gc_NULL(av);
    1278             :   /* smallest square root */
    1279       77748 :   av2 = avma; q = subii(p, a);
    1280       77748 :   if (cmpii(a, q) > 0) a = q; else set_avma(av2);
    1281       77750 :   return gerepileuptoint(av, a);
    1282             : }
    1283             : GEN
    1284     4834780 : Fp_sqrt(GEN a, GEN p) { return Fp_sqrt_i(a, NULL, p); }
    1285             : 
    1286             : /*********************************************************************/
    1287             : /**                        GCD & BEZOUT                             **/
    1288             : /*********************************************************************/
    1289             : 
    1290             : GEN
    1291    53379782 : lcmii(GEN x, GEN y)
    1292             : {
    1293             :   pari_sp av;
    1294             :   GEN a, b;
    1295    53379782 :   if (!signe(x) || !signe(y)) return gen_0;
    1296    53379797 :   av = avma; a = gcdii(x,y);
    1297    53378567 :   if (absequalii(a,y)) { set_avma(av); return absi(x); }
    1298    11691975 :   if (!equali1(a)) y = diviiexact(y,a);
    1299    11691967 :   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     5916507 : Fp_invgen(GEN x, GEN N, GEN *pd)
    1306             : {
    1307             :   GEN d, d0, e, v;
    1308     5916507 :   if (lgefint(N) == 3)
    1309             :   {
    1310     5131118 :     ulong dd, NN = N[2], xx = umodiu(x,NN);
    1311     5131134 :     if (!xx) { *pd = N; return gen_0; }
    1312     5131134 :     xx = Fl_invgen(xx, NN, &dd);
    1313     5132609 :     *pd = utoi(dd); return utoi(xx);
    1314             :   }
    1315      785389 :   *pd = d = bezout(x, N, &v, NULL);
    1316      785399 :   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      688507 :   e = diviiexact(N,d);
    1319      688507 :   d0 = Z_ppo(d, e); /* d = d0 d1, d0 coprime to N/d, rad(d1) | N/d */
    1320      688507 :   if (equali1(d0)) return v;
    1321      546142 :   if (!equalii(d,d0)) e = lcmii(e, diviiexact(d,d0));
    1322      546142 :   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     2415525 : gen_chinese(GEN x, GEN(*f)(GEN,GEN))
    1344             : {
    1345     2415525 :   GEN z = gassoc_proto(f,x,NULL);
    1346     2415517 :   if (z == gen_1) retmkintmod(gen_0,gen_1);
    1347     2415482 :   return z;
    1348             : }
    1349             : 
    1350             : GEN
    1351        2415 : chinese1(GEN x) { return gen_chinese(x,chinese); }
    1352             : 
    1353             : static GEN
    1354          21 : padic2mod(GEN x)
    1355             : {
    1356          21 :   pari_sp av = avma;
    1357          21 :   GEN pd = padic_pd(x), p = padic_p(x), u = padic_u(x);
    1358          21 :   long v = valp(x);
    1359          21 :   if (v < 0) pari_err_INV("chinese", mkintmod(gen_0, p));
    1360          21 :   if (v)
    1361             :   {
    1362           0 :     GEN pv = powiu(p, v);
    1363           0 :     pd = mulii(pd, pv);
    1364           0 :     u = mulii(u, pv);
    1365             :   }
    1366          21 :   return gerepilecopy(av, mkintmod(u, pd));
    1367             : 
    1368             : }
    1369             : /* x t_INTMOD, y t_POLMOD; promote x to t_POLMOD mod Pol(x.mod): makes Mod(0,1)
    1370             :  * a better "neutral" element */
    1371             : static GEN
    1372          21 : intmod2polmod(GEN x,GEN y)
    1373          21 : { retmkpolmod(gel(x,2), scalarpol_shallow(gel(x,1), varn(gel(y,1)))); }
    1374             : 
    1375             : GEN
    1376        5495 : chinese(GEN x, GEN y)
    1377             : {
    1378        5495 :   pari_sp av = avma;
    1379             :   long tx, ty;
    1380             :   GEN z;
    1381             : 
    1382        5495 :   if (!y) return chinese1(x);
    1383        5446 :   if (gidentical(x,y)) return gcopy(x);
    1384             :   /* allows GC optimization for this most frequent case */
    1385        5439 :   z = cgetg(3,t_INTMOD);
    1386        5439 :   tx = typ(x); if (tx == t_PADIC) { x = padic2mod(x); tx = t_INTMOD; }
    1387        5439 :   ty = typ(y); if (ty == t_PADIC) { y = padic2mod(y); ty = t_INTMOD; }
    1388        5439 :   if (tx == t_POLMOD && ty == t_INTMOD)
    1389          14 :   { y = intmod2polmod(y, x); ty = t_POLMOD; }
    1390        5439 :   if (ty == t_POLMOD && tx == t_INTMOD)
    1391           7 :   { x = intmod2polmod(x, y); tx = t_POLMOD; }
    1392        5439 :   if (tx == ty) switch(tx)
    1393             :   {
    1394        3892 :     case t_POLMOD:
    1395             :     {
    1396        3892 :       GEN A = gel(x,1), B = gel(y,1);
    1397        3892 :       GEN a = gel(x,2), b = gel(y,2), t, d, e, u, v;
    1398        3892 :       if (varn(A)!=varn(B)) pari_err_VAR("chinese",A,B);
    1399        3892 :       if (RgX_equal(A,B)) retmkpolmod(chinese(a,b), gcopy(A)); /*same modulus*/
    1400        3892 :       d = RgX_extgcd(A,B,&u,&v);
    1401        3892 :       e = gsub(b, a);
    1402        3892 :       if (!gequal0(gmod(e, d))) pari_err_OP("chinese",x,y);
    1403        3892 :       t = gdiv(A, d);
    1404        3892 :       e = gadd(a, gmul(gmul(u,t), e));
    1405             : 
    1406        3892 :       z = cgetg(3, t_POLMOD);
    1407        3892 :       gel(z,1) = RgX_mul(t, B);
    1408        3892 :       gel(z,2) = gmod(e, gel(z,1));
    1409        3892 :       return gerepileupto(av, z);
    1410             :     }
    1411        1519 :     case t_INTMOD:
    1412             :     {
    1413        1519 :       GEN A = gel(x,1), B = gel(y,1);
    1414        1519 :       GEN a = gel(x,2), b = gel(y,2), c, d, C, U;
    1415        1519 :       Z_chinese_pre(A, B, &C, &U, &d);
    1416        1519 :       c = Z_chinese_post(a, b, C, U, d);
    1417        1519 :       if (!c) pari_err_OP("chinese", x,y);
    1418        1519 :       set_avma((pari_sp)z); /* GC optimization */
    1419        1519 :       gel(z,1) = icopy(C);
    1420        1519 :       gel(z,2) = icopy(c); return z;
    1421             :     }
    1422          14 :     case t_POL:
    1423             :     {
    1424          14 :       long i, lx = lg(x), ly = lg(y);
    1425          14 :       if (varn(x) != varn(y)) pari_err_OP("chinese",x,y);
    1426          14 :       if (lx < ly) { swap(x,y); lswap(lx,ly); }
    1427          14 :       set_avma(av);
    1428          14 :       z = cgetg(lx, t_POL); z[1] = x[1];
    1429          42 :       for (i=2; i<ly; i++) gel(z,i) = chinese(gel(x,i),gel(y,i));
    1430          14 :       if (i < lx)
    1431             :       {
    1432          14 :         GEN _0 = Rg_get_0(y);
    1433          28 :         for (   ; i<lx; i++) gel(z,i) = chinese(gel(x,i),_0);
    1434             :       }
    1435          14 :       return z;
    1436             :     }
    1437          14 :     case t_VEC: case t_COL: case t_MAT:
    1438             :     {
    1439             :       long i, lx;
    1440          14 :       set_avma(av);
    1441          14 :       z = cgetg_copy(x, &lx); if (lx!=lg(y)) pari_err_OP("chinese",x,y);
    1442          42 :       for (i=1; i<lx; i++) gel(z,i) = chinese(gel(x,i),gel(y,i));
    1443          14 :       return z;
    1444             :     }
    1445             :   }
    1446           0 :   pari_err_OP("chinese",x,y);
    1447             :   return NULL; /* LCOV_EXCL_LINE */
    1448             : }
    1449             : 
    1450             : /* init chinese(Mod(.,A), Mod(.,B)) */
    1451             : void
    1452      271327 : Z_chinese_pre(GEN A, GEN B, GEN *pC, GEN *pU, GEN *pd)
    1453             : {
    1454      271327 :   GEN u, d = bezout(A,B,&u,NULL); /* U = u(A/d), u(A/d) + v(B/d) = 1 */
    1455      271329 :   GEN t = diviiexact(A,d);
    1456      271322 :   *pU = mulii(u, t);
    1457      271322 :   *pC = mulii(t, B); if (pd) *pd = d;
    1458      271322 : }
    1459             : /* Assume C = lcm(A, B), U = 0 mod (A/d), U = 1 mod (B/d), a = b mod d,
    1460             :  * where d = gcd(A,B) or NULL, return x = a (mod A), b (mod B).
    1461             :  * If d not NULL, check whether a = b mod d. */
    1462             : GEN
    1463     3010762 : Z_chinese_post(GEN a, GEN b, GEN C, GEN U, GEN d)
    1464             : {
    1465             :   GEN e;
    1466     3010762 :   if (!signe(a))
    1467             :   {
    1468      797085 :     if (d && !dvdii(b, d)) return NULL;
    1469      797085 :     return Fp_mul(b, U, C);
    1470             :   }
    1471     2213677 :   e = subii(b,a);
    1472     2213677 :   if (d && !dvdii(e, d)) return NULL;
    1473     2213677 :   return modii(addii(a, mulii(U, e)), C);
    1474             : }
    1475             : static ulong
    1476     1592118 : u_chinese_post(ulong a, ulong b, ulong C, ulong U)
    1477             : {
    1478     1592118 :   if (!a) return Fl_mul(b, U, C);
    1479     1592118 :   return Fl_add(a, Fl_mul(U, Fl_sub(b,a,C), C), C);
    1480             : }
    1481             : 
    1482             : GEN
    1483        2142 : Z_chinese(GEN a, GEN b, GEN A, GEN B)
    1484             : {
    1485        2142 :   pari_sp av = avma;
    1486        2142 :   GEN C, U; Z_chinese_pre(A, B, &C, &U, NULL);
    1487        2142 :   return gerepileuptoint(av, Z_chinese_post(a,b, C, U, NULL));
    1488             : }
    1489             : GEN
    1490      267609 : Z_chinese_all(GEN a, GEN b, GEN A, GEN B, GEN *pC)
    1491             : {
    1492      267609 :   GEN U; Z_chinese_pre(A, B, pC, &U, NULL);
    1493      267605 :   return Z_chinese_post(a,b, *pC, U, NULL);
    1494             : }
    1495             : 
    1496             : /* return lift(chinese(a mod A, b mod B))
    1497             :  * assume(A,B)=1, a,b,A,B integers and C = A*B */
    1498             : GEN
    1499      547401 : Z_chinese_coprime(GEN a, GEN b, GEN A, GEN B, GEN C)
    1500             : {
    1501      547401 :   pari_sp av = avma;
    1502      547401 :   GEN U = mulii(Fp_inv(A,B), A);
    1503      547402 :   return gerepileuptoint(av, Z_chinese_post(a,b,C,U, NULL));
    1504             : }
    1505             : ulong
    1506     1592112 : u_chinese_coprime(ulong a, ulong b, ulong A, ulong B, ulong C)
    1507     1592112 : { return u_chinese_post(a,b,C, A * Fl_inv(A % B,B)); }
    1508             : 
    1509             : /* chinese1 for coprime moduli in Z */
    1510             : static GEN
    1511     2191772 : chinese1_coprime_Z_aux(GEN x, GEN y)
    1512             : {
    1513     2191772 :   GEN z = cgetg(3, t_INTMOD);
    1514     2191772 :   GEN A = gel(x,1), a = gel(x, 2);
    1515     2191772 :   GEN B = gel(y,1), b = gel(y, 2), C = mulii(A,B);
    1516     2191772 :   pari_sp av = avma;
    1517     2191772 :   GEN U = mulii(Fp_inv(A,B), A);
    1518     2191772 :   gel(z,2) = gerepileuptoint(av, Z_chinese_post(a,b,C,U, NULL));
    1519     2191772 :   gel(z,1) = C; return z;
    1520             : }
    1521             : GEN
    1522     2413110 : chinese1_coprime_Z(GEN x) {return gen_chinese(x,chinese1_coprime_Z_aux);}
    1523             : 
    1524             : /*********************************************************************/
    1525             : /**                    MODULAR EXPONENTIATION                       **/
    1526             : /*********************************************************************/
    1527             : /* xa ZV or nv */
    1528             : GEN
    1529     2597408 : ZV_producttree(GEN xa)
    1530             : {
    1531     2597408 :   long n = lg(xa)-1;
    1532     2597408 :   long m = n==1 ? 1: expu(n-1)+1;
    1533     2597408 :   GEN T = cgetg(m+1, t_VEC), t;
    1534             :   long i, j, k;
    1535     2597407 :   t = cgetg(((n+1)>>1)+1, t_VEC);
    1536     2597402 :   if (typ(xa)==t_VECSMALL)
    1537             :   {
    1538     3472711 :     for (j=1, k=1; k<n; j++, k+=2)
    1539     2238563 :       gel(t, j) = muluu(xa[k], xa[k+1]);
    1540     1234148 :     if (k==n) gel(t, j) = utoi(xa[k]);
    1541             :   } else {
    1542     2823531 :     for (j=1, k=1; k<n; j++, k+=2)
    1543     1460269 :       gel(t, j) = mulii(gel(xa,k), gel(xa,k+1));
    1544     1363262 :     if (k==n) gel(t, j) = icopy(gel(xa,k));
    1545             :   }
    1546     2597410 :   gel(T,1) = t;
    1547     4145785 :   for (i=2; i<=m; i++)
    1548             :   {
    1549     1548375 :     GEN u = gel(T, i-1);
    1550     1548375 :     long n = lg(u)-1;
    1551     1548375 :     t = cgetg(((n+1)>>1)+1, t_VEC);
    1552     3481063 :     for (j=1, k=1; k<n; j++, k+=2)
    1553     1932688 :       gel(t, j) = mulii(gel(u, k), gel(u, k+1));
    1554     1548375 :     if (k==n) gel(t, j) = gel(u, k);
    1555     1548375 :     gel(T, i) = t;
    1556             :   }
    1557     2597410 :   return T;
    1558             : }
    1559             : 
    1560             : /* return [A mod P[i], i=1..#P], T = ZV_producttree(P) */
    1561             : GEN
    1562    57159375 : Z_ZV_mod_tree(GEN A, GEN P, GEN T)
    1563             : {
    1564             :   long i,j,k;
    1565    57159375 :   long m = lg(T)-1, n = lg(P)-1;
    1566             :   GEN t;
    1567    57159375 :   GEN Tp = cgetg(m+1, t_VEC);
    1568    57109788 :   gel(Tp, m) = mkvec(modii(A, gmael(T,m,1)));
    1569   119130154 :   for (i=m-1; i>=1; i--)
    1570             :   {
    1571    62101185 :     GEN u = gel(T, i);
    1572    62101185 :     GEN v = gel(Tp, i+1);
    1573    62101185 :     long n = lg(u)-1;
    1574    62101185 :     t = cgetg(n+1, t_VEC);
    1575   148494517 :     for (j=1, k=1; k<n; j++, k+=2)
    1576             :     {
    1577    86471456 :       gel(t, k)   = modii(gel(v, j), gel(u, k));
    1578    86531720 :       gel(t, k+1) = modii(gel(v, j), gel(u, k+1));
    1579             :     }
    1580    62023061 :     if (k==n) gel(t, k) = gel(v, j);
    1581    62023061 :     gel(Tp, i) = t;
    1582             :   }
    1583             :   {
    1584    57028969 :     GEN u = gel(T, i+1);
    1585    57028969 :     GEN v = gel(Tp, i+1);
    1586    57028969 :     long l = lg(u)-1;
    1587    57028969 :     if (typ(P)==t_VECSMALL)
    1588             :     {
    1589    54434676 :       GEN R = cgetg(n+1, t_VECSMALL);
    1590   194655159 :       for (j=1, k=1; j<=l; j++, k+=2)
    1591             :       {
    1592   139944575 :         uel(R,k) = umodiu(gel(v, j), P[k]);
    1593   140202512 :         if (k < n)
    1594   110691753 :           uel(R,k+1) = umodiu(gel(v, j), P[k+1]);
    1595             :       }
    1596    54710584 :       return R;
    1597             :     }
    1598             :     else
    1599             :     {
    1600     2594293 :       GEN R = cgetg(n+1, t_VEC);
    1601     7123807 :       for (j=1, k=1; j<=l; j++, k+=2)
    1602             :       {
    1603     4526758 :         gel(R,k) = modii(gel(v, j), gel(P,k));
    1604     4526764 :         if (k < n)
    1605     3695611 :           gel(R,k+1) = modii(gel(v, j), gel(P,k+1));
    1606             :       }
    1607     2597049 :       return R;
    1608             :     }
    1609             :   }
    1610             : }
    1611             : 
    1612             : /* T = ZV_producttree(P), R = ZV_chinesetree(P,T) */
    1613             : GEN
    1614    39649264 : ZV_chinese_tree(GEN A, GEN P, GEN T, GEN R)
    1615             : {
    1616    39649264 :   long m = lg(T)-1, n = lg(A)-1;
    1617             :   long i,j,k;
    1618    39649264 :   GEN Tp = cgetg(m+1, t_VEC);
    1619    39641762 :   GEN M = gel(T, 1);
    1620    39641762 :   GEN t = cgetg(lg(M), t_VEC);
    1621    39587837 :   if (typ(P)==t_VECSMALL)
    1622             :   {
    1623    83870007 :     for (j=1, k=1; k<n; j++, k+=2)
    1624             :     {
    1625    60872294 :       pari_sp av = avma;
    1626    60872294 :       GEN a = mului(A[k], gel(R,k)), b = mului(A[k+1], gel(R,k+1));
    1627    60762806 :       GEN tj = modii(addii(mului(P[k],b), mului(P[k+1],a)), gel(M,j));
    1628    60870466 :       gel(t, j) = gerepileuptoint(av, tj);
    1629             :     }
    1630    22997713 :     if (k==n) gel(t, j) = modii(mului(A[k], gel(R,k)), gel(M, j));
    1631             :   } else
    1632             :   {
    1633    35264578 :     for (j=1, k=1; k<n; j++, k+=2)
    1634             :     {
    1635    18638551 :       pari_sp av = avma;
    1636    18638551 :       GEN a = mulii(gel(A,k), gel(R,k)), b = mulii(gel(A,k+1), gel(R,k+1));
    1637    18639863 :       GEN tj = modii(addii(mulii(gel(P,k),b), mulii(gel(P,k+1),a)), gel(M,j));
    1638    18687247 :       gel(t, j) = gerepileuptoint(av, tj);
    1639             :     }
    1640    16626027 :     if (k==n) gel(t, j) = modii(mulii(gel(A,k), gel(R,k)), gel(M, j));
    1641             :   }
    1642    39618062 :   gel(Tp, 1) = t;
    1643    73922187 :   for (i=2; i<=m; i++)
    1644             :   {
    1645    34280559 :     GEN u = gel(T, i-1), M = gel(T, i);
    1646    34280559 :     GEN t = cgetg(lg(M), t_VEC);
    1647    34278669 :     GEN v = gel(Tp, i-1);
    1648    34278669 :     long n = lg(v)-1;
    1649    89849272 :     for (j=1, k=1; k<n; j++, k+=2)
    1650             :     {
    1651    55545147 :       pari_sp av = avma;
    1652    55510708 :       gel(t, j) = gerepileuptoint(av, modii(addii(mulii(gel(u, k), gel(v, k+1)),
    1653    55545147 :             mulii(gel(u, k+1), gel(v, k))), gel(M, j)));
    1654             :     }
    1655    34304125 :     if (k==n) gel(t, j) = gel(v, k);
    1656    34304125 :     gel(Tp, i) = t;
    1657             :   }
    1658    39641628 :   return gmael(Tp,m,1);
    1659             : }
    1660             : 
    1661             : static GEN
    1662     1525268 : ncV_polint_center_tree(GEN vA, GEN P, GEN T, GEN R, GEN m2)
    1663             : {
    1664     1525268 :   long i, l = lg(gel(vA,1)), n = lg(P);
    1665     1525268 :   GEN mod = gmael(T, lg(T)-1, 1), V = cgetg(l, t_COL);
    1666    33856792 :   for (i=1; i < l; i++)
    1667             :   {
    1668    32331911 :     pari_sp av = avma;
    1669    32331911 :     GEN c, A = cgetg(n, typ(P));
    1670             :     long j;
    1671   186127849 :     for (j=1; j < n; j++) A[j] = mael(vA,j,i);
    1672    32297561 :     c = Fp_center(ZV_chinese_tree(A, P, T, R), mod, m2);
    1673    32334497 :     gel(V,i) = gerepileuptoint(av, c);
    1674             :   }
    1675     1524881 :   return V;
    1676             : }
    1677             : 
    1678             : static GEN
    1679      722673 : nxV_polint_center_tree(GEN vA, GEN P, GEN T, GEN R, GEN m2)
    1680             : {
    1681      722673 :   long i, j, l, n = lg(P);
    1682      722673 :   GEN mod = gmael(T, lg(T)-1, 1), V, w;
    1683      722673 :   w = cgetg(n, t_VECSMALL);
    1684     2552923 :   for(j=1; j<n; j++) w[j] = lg(gel(vA,j));
    1685      722672 :   l = vecsmall_max(w);
    1686      722670 :   V = cgetg(l, t_POL);
    1687      722638 :   V[1] = mael(vA,1,1);
    1688     5310589 :   for (i=2; i < l; i++)
    1689             :   {
    1690     4587927 :     pari_sp av = avma;
    1691     4587927 :     GEN c, A = cgetg(n, typ(P));
    1692     4587461 :     if (typ(P)==t_VECSMALL)
    1693    12176281 :       for (j=1; j < n; j++) A[j] = i < w[j] ? mael(vA,j,i): 0;
    1694             :     else
    1695     5791603 :       for (j=1; j < n; j++) gel(A,j) = i < w[j] ? gmael(vA,j,i): gen_0;
    1696     4587461 :     c = Fp_center(ZV_chinese_tree(A, P, T, R), mod, m2);
    1697     4588280 :     gel(V,i) = gerepileuptoint(av, c);
    1698             :   }
    1699      722662 :   return ZX_renormalize(V, l);
    1700             : }
    1701             : 
    1702             : static GEN
    1703        4620 : nxCV_polint_center_tree(GEN vA, GEN P, GEN T, GEN R, GEN m2)
    1704             : {
    1705        4620 :   long i, j, l = lg(gel(vA,1)), n = lg(P);
    1706        4620 :   GEN A = cgetg(n, t_VEC);
    1707        4620 :   GEN V = cgetg(l, t_COL);
    1708       90958 :   for (i=1; i < l; i++)
    1709             :   {
    1710      335252 :     for (j=1; j < n; j++) gel(A,j) = gmael(vA,j,i);
    1711       86338 :     gel(V,i) = nxV_polint_center_tree(A, P, T, R, m2);
    1712             :   }
    1713        4620 :   return V;
    1714             : }
    1715             : 
    1716             : static GEN
    1717      419982 : polint_chinese(GEN worker, GEN mA, GEN P)
    1718             : {
    1719      419982 :   long cnt, pending, n, i, j, l = lg(gel(mA,1));
    1720             :   struct pari_mt pt;
    1721             :   GEN done, va, M, A;
    1722             :   pari_timer ti;
    1723             : 
    1724      419982 :   if (l == 1) return cgetg(1, t_MAT);
    1725      390930 :   cnt = pending = 0;
    1726      390930 :   n = lg(P);
    1727      390930 :   A = cgetg(n, t_VEC);
    1728      390930 :   va = mkvec(A);
    1729      390930 :   M = cgetg(l, t_MAT);
    1730      390930 :   if (DEBUGLEVEL>4) timer_start(&ti);
    1731      390930 :   if (DEBUGLEVEL>5) err_printf("Start parallel Chinese remainder: ");
    1732      390930 :   mt_queue_start_lim(&pt, worker, l-1);
    1733     1420482 :   for (i=1; i<l || pending; i++)
    1734             :   {
    1735             :     long workid;
    1736     3863528 :     for(j=1; j < n; j++) gel(A,j) = gmael(mA,j,i);
    1737     1029552 :     mt_queue_submit(&pt, i, i<l? va: NULL);
    1738     1029552 :     done = mt_queue_get(&pt, &workid, &pending);
    1739     1029552 :     if (done)
    1740             :     {
    1741      989802 :       gel(M,workid) = done;
    1742      989802 :       if (DEBUGLEVEL>5) err_printf("%ld%% ",(++cnt)*100/(l-1));
    1743             :     }
    1744             :   }
    1745      390930 :   if (DEBUGLEVEL>5) err_printf("\n");
    1746      390930 :   if (DEBUGLEVEL>4) timer_printf(&ti, "nmV_chinese_center");
    1747      390930 :   mt_queue_end(&pt);
    1748      390930 :   return M;
    1749             : }
    1750             : 
    1751             : GEN
    1752         840 : nxMV_polint_center_tree_worker(GEN vA, GEN T, GEN R, GEN P, GEN m2)
    1753             : {
    1754         840 :   return nxCV_polint_center_tree(vA, P, T, R, m2);
    1755             : }
    1756             : 
    1757             : static GEN
    1758         431 : nxMV_polint_center_tree_seq(GEN vA, GEN P, GEN T, GEN R, GEN m2)
    1759             : {
    1760         431 :   long i, j, l = lg(gel(vA,1)), n = lg(P);
    1761         431 :   GEN A = cgetg(n, t_VEC);
    1762         431 :   GEN V = cgetg(l, t_MAT);
    1763        4211 :   for (i=1; i < l; i++)
    1764             :   {
    1765       15317 :     for (j=1; j < n; j++) gel(A,j) = gmael(vA,j,i);
    1766        3780 :     gel(V,i) = nxCV_polint_center_tree(A, P, T, R, m2);
    1767             :   }
    1768         431 :   return V;
    1769             : }
    1770             : 
    1771             : static GEN
    1772          90 : nxMV_polint_center_tree(GEN mA, GEN P, GEN T, GEN R, GEN m2)
    1773             : {
    1774          90 :   GEN worker = snm_closure(is_entry("_nxMV_polint_worker"), mkvec4(T, R, P, m2));
    1775          90 :   return polint_chinese(worker, mA, P);
    1776             : }
    1777             : 
    1778             : static GEN
    1779      141704 : nmV_polint_center_tree_seq(GEN vA, GEN P, GEN T, GEN R, GEN m2)
    1780             : {
    1781      141704 :   long i, j, l = lg(gel(vA,1)), n = lg(P);
    1782      141704 :   GEN A = cgetg(n, t_VEC);
    1783      141704 :   GEN V = cgetg(l, t_MAT);
    1784      662311 :   for (i=1; i < l; i++)
    1785             :   {
    1786     2985491 :     for (j=1; j < n; j++) gel(A,j) = gmael(vA,j,i);
    1787      520608 :     gel(V,i) = ncV_polint_center_tree(A, P, T, R, m2);
    1788             :   }
    1789      141703 :   return V;
    1790             : }
    1791             : 
    1792             : GEN
    1793      988931 : nmV_polint_center_tree_worker(GEN vA, GEN T, GEN R, GEN P, GEN m2)
    1794             : {
    1795      988931 :   return ncV_polint_center_tree(vA, P, T, R, m2);
    1796             : }
    1797             : 
    1798             : static GEN
    1799      419892 : nmV_polint_center_tree(GEN mA, GEN P, GEN T, GEN R, GEN m2)
    1800             : {
    1801      419892 :   GEN worker = snm_closure(is_entry("_polint_worker"), mkvec4(T, R, P, m2));
    1802      419892 :   return polint_chinese(worker, mA, P);
    1803             : }
    1804             : 
    1805             : /* return [A mod P[i], i=1..#P] */
    1806             : GEN
    1807           0 : Z_ZV_mod(GEN A, GEN P)
    1808             : {
    1809           0 :   pari_sp av = avma;
    1810           0 :   return gerepilecopy(av, Z_ZV_mod_tree(A, P, ZV_producttree(P)));
    1811             : }
    1812             : /* P a t_VECSMALL */
    1813             : GEN
    1814           0 : Z_nv_mod(GEN A, GEN P)
    1815             : {
    1816           0 :   pari_sp av = avma;
    1817           0 :   return gerepileuptoleaf(av, Z_ZV_mod_tree(A, P, ZV_producttree(P)));
    1818             : }
    1819             : /* B a ZX, T = ZV_producttree(P) */
    1820             : GEN
    1821     2401293 : ZX_nv_mod_tree(GEN B, GEN A, GEN T)
    1822             : {
    1823             :   pari_sp av;
    1824     2401293 :   long i, j, l = lg(B), n = lg(A)-1;
    1825     2401293 :   GEN V = cgetg(n+1, t_VEC);
    1826    11404198 :   for (j=1; j <= n; j++)
    1827             :   {
    1828     9003083 :     gel(V, j) = cgetg(l, t_VECSMALL);
    1829     9002940 :     mael(V, j, 1) = B[1]&VARNBITS;
    1830             :   }
    1831     2401115 :   av = avma;
    1832    15677920 :   for (i=2; i < l; i++)
    1833             :   {
    1834    13278313 :     GEN v = Z_ZV_mod_tree(gel(B, i), A, T);
    1835    87814271 :     for (j=1; j <= n; j++)
    1836    74546543 :       mael(V, j, i) = v[j];
    1837    13267728 :     set_avma(av);
    1838             :   }
    1839    11402916 :   for (j=1; j <= n; j++)
    1840     9003273 :     (void) Flx_renormalize(gel(V, j), l);
    1841     2399643 :   return V;
    1842             : }
    1843             : 
    1844             : static GEN
    1845     1191427 : to_ZX(GEN a, long v) { return typ(a)==t_INT? scalarpol(a,v): a; }
    1846             : 
    1847             : GEN
    1848       86827 : ZXX_nv_mod_tree(GEN P, GEN xa, GEN T, long w)
    1849             : {
    1850       86827 :   pari_sp av = avma;
    1851       86827 :   long i, j, l = lg(P), n = lg(xa)-1;
    1852       86827 :   GEN V = cgetg(n+1, t_VEC);
    1853      374755 :   for (j=1; j <= n; j++)
    1854             :   {
    1855      287928 :     gel(V, j) = cgetg(l, t_POL);
    1856      287928 :     mael(V, j, 1) = P[1]&VARNBITS;
    1857             :   }
    1858     1197022 :   for (i=2; i < l; i++)
    1859             :   {
    1860     1110196 :     GEN v = ZX_nv_mod_tree(to_ZX(gel(P, i), w), xa, T);
    1861     4841537 :     for (j=1; j <= n; j++)
    1862     3731342 :       gmael(V, j, i) = gel(v,j);
    1863             :   }
    1864      374754 :   for (j=1; j <= n; j++)
    1865      287928 :     (void) FlxX_renormalize(gel(V, j), l);
    1866       86826 :   return gerepilecopy(av, V);
    1867             : }
    1868             : 
    1869             : GEN
    1870        4054 : ZXC_nv_mod_tree(GEN C, GEN xa, GEN T, long w)
    1871             : {
    1872        4054 :   pari_sp av = avma;
    1873        4054 :   long i, j, l = lg(C), n = lg(xa)-1;
    1874        4054 :   GEN V = cgetg(n+1, t_VEC);
    1875       16970 :   for (j = 1; j <= n; j++)
    1876       12916 :     gel(V, j) = cgetg(l, t_COL);
    1877       85280 :   for (i = 1; i < l; i++)
    1878             :   {
    1879       81228 :     GEN v = ZX_nv_mod_tree(to_ZX(gel(C, i), w), xa, T);
    1880      359623 :     for (j = 1; j <= n; j++)
    1881      278397 :       gmael(V, j, i) = gel(v,j);
    1882             :   }
    1883        4052 :   return gerepilecopy(av, V);
    1884             : }
    1885             : 
    1886             : GEN
    1887         431 : ZXM_nv_mod_tree(GEN M, GEN xa, GEN T, long w)
    1888             : {
    1889         431 :   pari_sp av = avma;
    1890         431 :   long i, j, l = lg(M), n = lg(xa)-1;
    1891         431 :   GEN V = cgetg(n+1, t_VEC);
    1892        2086 :   for (j=1; j <= n; j++)
    1893        1655 :     gel(V, j) = cgetg(l, t_MAT);
    1894        4211 :   for (i=1; i < l; i++)
    1895             :   {
    1896        3780 :     GEN v = ZXC_nv_mod_tree(gel(M, i), xa, T, w);
    1897       15317 :     for (j=1; j <= n; j++)
    1898       11537 :       gmael(V, j, i) = gel(v,j);
    1899             :   }
    1900         431 :   return gerepilecopy(av, V);
    1901             : }
    1902             : 
    1903             : GEN
    1904     1273764 : ZV_nv_mod_tree(GEN B, GEN A, GEN T)
    1905             : {
    1906             :   pari_sp av;
    1907     1273764 :   long i, j, l = lg(B), n = lg(A)-1;
    1908     1273764 :   GEN V = cgetg(n+1, t_VEC);
    1909     6533905 :   for (j=1; j <= n; j++) gel(V, j) = cgetg(l, t_VECSMALL);
    1910     1273669 :   av = avma;
    1911    42488870 :   for (i=1; i < l; i++)
    1912             :   {
    1913    41222240 :     GEN v = Z_ZV_mod_tree(gel(B, i), A, T);
    1914   218849642 :     for (j=1; j <= n; j++) mael(V, j, i) = v[j];
    1915    41162381 :     set_avma(av);
    1916             :   }
    1917     1266630 :   return V;
    1918             : }
    1919             : 
    1920             : static GEN
    1921      241416 : ZM_nv_mod_tree_t(GEN M, GEN xa, GEN T, long t)
    1922             : {
    1923      241416 :   pari_sp av = avma;
    1924      241416 :   long i, j, l = lg(M), n = lg(xa)-1;
    1925      241416 :   GEN V = cgetg(n+1, t_VEC);
    1926     1322619 :   for (j=1; j <= n; j++) gel(V, j) = cgetg(l, t);
    1927     1514951 :   for (i=1; i < l; i++)
    1928             :   {
    1929     1273538 :     GEN v = ZV_nv_mod_tree(gel(M, i), xa, T);
    1930     6533427 :     for (j=1; j <= n; j++) gmael(V, j, i) = gel(v,j);
    1931             :   }
    1932      241413 :   return gerepilecopy(av, V);
    1933             : }
    1934             : 
    1935             : GEN
    1936      235959 : ZM_nv_mod_tree(GEN M, GEN xa, GEN T)
    1937      235959 : { return ZM_nv_mod_tree_t(M, xa, T, t_MAT); }
    1938             : 
    1939             : GEN
    1940        5457 : ZVV_nv_mod_tree(GEN M, GEN xa, GEN T)
    1941        5457 : { return ZM_nv_mod_tree_t(M, xa, T, t_VEC); }
    1942             : 
    1943             : static GEN
    1944     2593199 : ZV_sqr(GEN z)
    1945             : {
    1946     2593199 :   long i,l = lg(z);
    1947     2593199 :   GEN x = cgetg(l, t_VEC);
    1948     2593196 :   if (typ(z)==t_VECSMALL)
    1949     6163227 :     for (i=1; i<l; i++) gel(x,i) = sqru(z[i]);
    1950             :   else
    1951     4626615 :     for (i=1; i<l; i++) gel(x,i) = sqri(gel(z,i));
    1952     2593175 :   return x;
    1953             : }
    1954             : 
    1955             : static GEN
    1956    13412971 : ZT_sqr(GEN x)
    1957             : {
    1958    13412971 :   if (typ(x) == t_INT) return sqri(x);
    1959    17548072 :   pari_APPLY_type(t_VEC, ZT_sqr(gel(x,i)))
    1960             : }
    1961             : 
    1962             : static GEN
    1963     2593194 : ZV_invdivexact(GEN y, GEN x)
    1964             : {
    1965     2593194 :   long i, l = lg(y);
    1966     2593194 :   GEN z = cgetg(l,t_VEC);
    1967     2593190 :   if (typ(x)==t_VECSMALL)
    1968     6163091 :     for (i=1; i<l; i++)
    1969             :     {
    1970     4929302 :       pari_sp av = avma;
    1971     4929302 :       ulong a = Fl_inv(umodiu(diviuexact(gel(y,i),x[i]), x[i]), x[i]);
    1972     4929506 :       set_avma(av); gel(z,i) = utoi(a);
    1973             :     }
    1974             :   else
    1975     4626618 :     for (i=1; i<l; i++)
    1976     3267217 :       gel(z,i) = Fp_inv(diviiexact(gel(y,i), gel(x,i)), gel(x,i));
    1977     2593190 :   return z;
    1978             : }
    1979             : 
    1980             : /* P t_VECSMALL or t_VEC of t_INT  */
    1981             : GEN
    1982     2593187 : ZV_chinesetree(GEN P, GEN T)
    1983             : {
    1984     2593187 :   GEN T2 = ZT_sqr(T), P2 = ZV_sqr(P);
    1985     2593188 :   GEN mod = gmael(T,lg(T)-1,1);
    1986     2593188 :   return ZV_invdivexact(Z_ZV_mod_tree(mod, P2, T2), P);
    1987             : }
    1988             : 
    1989             : static GEN
    1990     1010104 : gc_chinese(pari_sp av, GEN T, GEN a, GEN *pt_mod)
    1991             : {
    1992     1010104 :   if (!pt_mod)
    1993       12377 :     return gerepileupto(av, a);
    1994             :   else
    1995             :   {
    1996      997727 :     GEN mod = gmael(T, lg(T)-1, 1);
    1997      997727 :     gerepileall(av, 2, &a, &mod);
    1998      997727 :     *pt_mod = mod;
    1999      997727 :     return a;
    2000             :   }
    2001             : }
    2002             : 
    2003             : GEN
    2004      157264 : ZV_chinese_center(GEN A, GEN P, GEN *pt_mod)
    2005             : {
    2006      157264 :   pari_sp av = avma;
    2007      157264 :   GEN T = ZV_producttree(P);
    2008      157264 :   GEN R = ZV_chinesetree(P, T);
    2009      157264 :   GEN a = ZV_chinese_tree(A, P, T, R);
    2010      157264 :   GEN mod = gmael(T, lg(T)-1, 1);
    2011      157264 :   GEN ca = Fp_center(a, mod, shifti(mod,-1));
    2012      157264 :   return gc_chinese(av, T, ca, pt_mod);
    2013             : }
    2014             : 
    2015             : GEN
    2016        5141 : ZV_chinese(GEN A, GEN P, GEN *pt_mod)
    2017             : {
    2018        5141 :   pari_sp av = avma;
    2019        5141 :   GEN T = ZV_producttree(P);
    2020        5141 :   GEN R = ZV_chinesetree(P, T);
    2021        5141 :   GEN a = ZV_chinese_tree(A, P, T, R);
    2022        5141 :   return gc_chinese(av, T, a, pt_mod);
    2023             : }
    2024             : 
    2025             : GEN
    2026      218876 : nxV_chinese_center_tree(GEN A, GEN P, GEN T, GEN R)
    2027             : {
    2028      218876 :   pari_sp av = avma;
    2029      218876 :   GEN m2 = shifti(gmael(T, lg(T)-1, 1), -1);
    2030      218875 :   GEN a = nxV_polint_center_tree(A, P, T, R, m2);
    2031      218876 :   return gerepileupto(av, a);
    2032             : }
    2033             : 
    2034             : GEN
    2035      417453 : nxV_chinese_center(GEN A, GEN P, GEN *pt_mod)
    2036             : {
    2037      417453 :   pari_sp av = avma;
    2038      417453 :   GEN T = ZV_producttree(P);
    2039      417454 :   GEN R = ZV_chinesetree(P, T);
    2040      417452 :   GEN m2 = shifti(gmael(T, lg(T)-1, 1), -1);
    2041      417452 :   GEN a = nxV_polint_center_tree(A, P, T, R, m2);
    2042      417454 :   return gc_chinese(av, T, a, pt_mod);
    2043             : }
    2044             : 
    2045             : GEN
    2046       10263 : ncV_chinese_center(GEN A, GEN P, GEN *pt_mod)
    2047             : {
    2048       10263 :   pari_sp av = avma;
    2049       10263 :   GEN T = ZV_producttree(P);
    2050       10263 :   GEN R = ZV_chinesetree(P, T);
    2051       10263 :   GEN m2 = shifti(gmael(T, lg(T)-1, 1), -1);
    2052       10263 :   GEN a = ncV_polint_center_tree(A, P, T, R, m2);
    2053       10263 :   return gc_chinese(av, T, a, pt_mod);
    2054             : }
    2055             : 
    2056             : GEN
    2057        5457 : ncV_chinese_center_tree(GEN A, GEN P, GEN T, GEN R)
    2058             : {
    2059        5457 :   pari_sp av = avma;
    2060        5457 :   GEN m2 = shifti(gmael(T, lg(T)-1, 1), -1);
    2061        5457 :   GEN a = ncV_polint_center_tree(A, P, T, R, m2);
    2062        5457 :   return gerepileupto(av, a);
    2063             : }
    2064             : 
    2065             : GEN
    2066           0 : nmV_chinese_center_tree(GEN A, GEN P, GEN T, GEN R)
    2067             : {
    2068           0 :   pari_sp av = avma;
    2069           0 :   GEN m2 = shifti(gmael(T, lg(T)-1, 1), -1);
    2070           0 :   GEN a = nmV_polint_center_tree(A, P, T, R, m2);
    2071           0 :   return gerepileupto(av, a);
    2072             : }
    2073             : 
    2074             : GEN
    2075      141704 : nmV_chinese_center_tree_seq(GEN A, GEN P, GEN T, GEN R)
    2076             : {
    2077      141704 :   pari_sp av = avma;
    2078      141704 :   GEN m2 = shifti(gmael(T, lg(T)-1, 1), -1);
    2079      141704 :   GEN a = nmV_polint_center_tree_seq(A, P, T, R, m2);
    2080      141703 :   return gerepileupto(av, a);
    2081             : }
    2082             : 
    2083             : GEN
    2084      419892 : nmV_chinese_center(GEN A, GEN P, GEN *pt_mod)
    2085             : {
    2086      419892 :   pari_sp av = avma;
    2087      419892 :   GEN T = ZV_producttree(P);
    2088      419892 :   GEN R = ZV_chinesetree(P, T);
    2089      419892 :   GEN m2 = shifti(gmael(T, lg(T)-1, 1), -1);
    2090      419892 :   GEN a = nmV_polint_center_tree(A, P, T, R, m2);
    2091      419892 :   return gc_chinese(av, T, a, pt_mod);
    2092             : }
    2093             : 
    2094             : GEN
    2095           0 : nxCV_chinese_center_tree(GEN A, GEN P, GEN T, GEN R)
    2096             : {
    2097           0 :   pari_sp av = avma;
    2098           0 :   GEN m2 = shifti(gmael(T, lg(T)-1, 1), -1);
    2099           0 :   GEN a = nxCV_polint_center_tree(A, P, T, R, m2);
    2100           0 :   return gerepileupto(av, a);
    2101             : }
    2102             : 
    2103             : GEN
    2104           0 : nxCV_chinese_center(GEN A, GEN P, GEN *pt_mod)
    2105             : {
    2106           0 :   pari_sp av = avma;
    2107           0 :   GEN T = ZV_producttree(P);
    2108           0 :   GEN R = ZV_chinesetree(P, T);
    2109           0 :   GEN m2 = shifti(gmael(T, lg(T)-1, 1), -1);
    2110           0 :   GEN a = nxCV_polint_center_tree(A, P, T, R, m2);
    2111           0 :   return gc_chinese(av, T, a, pt_mod);
    2112             : }
    2113             : 
    2114             : GEN
    2115         431 : nxMV_chinese_center_tree_seq(GEN A, GEN P, GEN T, GEN R)
    2116             : {
    2117         431 :   pari_sp av = avma;
    2118         431 :   GEN m2 = shifti(gmael(T, lg(T)-1, 1), -1);
    2119         431 :   GEN a = nxMV_polint_center_tree_seq(A, P, T, R, m2);
    2120         431 :   return gerepileupto(av, a);
    2121             : }
    2122             : 
    2123             : GEN
    2124          90 : nxMV_chinese_center(GEN A, GEN P, GEN *pt_mod)
    2125             : {
    2126          90 :   pari_sp av = avma;
    2127          90 :   GEN T = ZV_producttree(P);
    2128          90 :   GEN R = ZV_chinesetree(P, T);
    2129          90 :   GEN m2 = shifti(gmael(T, lg(T)-1, 1), -1);
    2130          90 :   GEN a = nxMV_polint_center_tree(A, P, T, R, m2);
    2131          90 :   return gc_chinese(av, T, a, pt_mod);
    2132             : }
    2133             : 
    2134             : /**********************************************************************
    2135             :  **                    Powering  over (Z/NZ)^*, small N              **
    2136             :  **********************************************************************/
    2137             : 
    2138             : /* 2^n mod p; assume n > 1 */
    2139             : static ulong
    2140    12488581 : Fl_2powu_pre(ulong n, ulong p, ulong pi)
    2141             : {
    2142    12488581 :   ulong y = 2;
    2143    12488581 :   int j = 1+bfffo(n);
    2144             :   /* normalize, i.e set highest bit to 1 (we know n != 0) */
    2145    12488581 :   n<<=j; j = BITS_IN_LONG-j; /* first bit is now implicit */
    2146   559928226 :   for (; j; n<<=1,j--)
    2147             :   {
    2148   547459044 :     y = Fl_sqr_pre(y,p,pi);
    2149   547457367 :     if (n & HIGHBIT) y = Fl_double(y, p);
    2150             :   }
    2151    12469182 :   return y;
    2152             : }
    2153             : 
    2154             : /* 2^n mod p; assume n > 1 and !(p & HIGHMASK) */
    2155             : static ulong
    2156     4308771 : Fl_2powu(ulong n, ulong p)
    2157             : {
    2158     4308771 :   ulong y = 2;
    2159     4308771 :   int j = 1+bfffo(n);
    2160             :   /* normalize, i.e set highest bit to 1 (we know n != 0) */
    2161     4308771 :   n<<=j; j = BITS_IN_LONG-j; /* first bit is now implicit */
    2162    25453726 :   for (; j; n<<=1,j--)
    2163             :   {
    2164    21144910 :     y = (y*y) % p;
    2165    21144910 :     if (n & HIGHBIT) y = Fl_double(y, p);
    2166             :   }
    2167     4308816 :   return y;
    2168             : }
    2169             : 
    2170             : /* allow pi = 0 */
    2171             : ulong
    2172   151583734 : Fl_powu_pre(ulong x, ulong n0, ulong p, ulong pi)
    2173             : {
    2174             :   ulong y, z, n;
    2175   151583734 :   if (!pi) return Fl_powu(x, n0, p);
    2176   149134695 :   if (n0 <= 1)
    2177             :   { /* frequent special cases */
    2178    10182791 :     if (n0 == 1) return x;
    2179      105215 :     if (n0 == 0) return 1;
    2180             :   }
    2181   138951900 :   if (x <= 2)
    2182             :   {
    2183    12748019 :     if (x == 2) return Fl_2powu_pre(n0, p, pi);
    2184      258624 :     return x; /* 0 or 1 */
    2185             :   }
    2186   126203881 :   y = 1; z = x; n = n0;
    2187             :   for(;;)
    2188             :   {
    2189   644960094 :     if (n&1) y = Fl_mul_pre(y,z,p,pi);
    2190   645311959 :     n>>=1; if (!n) return y;
    2191   518681587 :     z = Fl_sqr_pre(z,p,pi);
    2192             :   }
    2193             : }
    2194             : 
    2195             : ulong
    2196   139877929 : Fl_powu(ulong x, ulong n0, ulong p)
    2197             : {
    2198             :   ulong y, z, n;
    2199   139877929 :   if (n0 <= 2)
    2200             :   { /* frequent special cases */
    2201    66032090 :     if (n0 == 2) return Fl_sqr(x,p);
    2202    32407571 :     if (n0 == 1) return x;
    2203     1968366 :     if (n0 == 0) return 1;
    2204             :   }
    2205    73817077 :   if (x <= 1) return x; /* 0 or 1 */
    2206    73378692 :   if (p & HIGHMASK)
    2207     7903633 :     return Fl_powu_pre(x, n0, p, get_Fl_red(p));
    2208    65475059 :   if (x == 2) return Fl_2powu(n0, p);
    2209    61166293 :   y = 1; z = x; n = n0;
    2210             :   for(;;)
    2211             :   {
    2212   262042527 :     if (n&1) y = (y*z) % p;
    2213   262042527 :     n>>=1; if (!n) return y;
    2214   200876234 :     z = (z*z) % p;
    2215             :   }
    2216             : }
    2217             : 
    2218             : /* Reduce data dependency to maximize internal parallelism; allow pi = 0 */
    2219             : GEN
    2220    12804605 : Fl_powers_pre(ulong x, long n, ulong p, ulong pi)
    2221             : {
    2222             :   long i, k;
    2223    12804605 :   GEN z = cgetg(n + 2, t_VECSMALL);
    2224    12798066 :   z[1] = 1; if (n == 0) return z;
    2225    12798066 :   z[2] = x;
    2226    12798066 :   if (pi)
    2227             :   {
    2228    90107016 :     for (i = 3, k=2; i <= n; i+=2, k++)
    2229             :     {
    2230    77517612 :       z[i] = Fl_sqr_pre(z[k], p, pi);
    2231    77527406 :       z[i+1] = Fl_mul_pre(z[k], z[k+1], p, pi);
    2232             :     }
    2233    12589404 :     if (i==n+1) z[i] = Fl_sqr_pre(z[k], p, pi);
    2234             :   }
    2235      213054 :   else if (p & HIGHMASK)
    2236             :   {
    2237           0 :     for (i = 3, k=2; i <= n; i+=2, k++)
    2238             :     {
    2239           0 :       z[i] = Fl_sqr(z[k], p);
    2240           0 :       z[i+1] = Fl_mul(z[k], z[k+1], p);
    2241             :     }
    2242           0 :     if (i==n+1) z[i] = Fl_sqr(z[k], p);
    2243             :   }
    2244             :   else
    2245   400545629 :     for (i = 2; i <= n; i++) z[i+1] = (z[i] * x) % p;
    2246    12804631 :   return z;
    2247             : }
    2248             : 
    2249             : GEN
    2250      296078 : Fl_powers(ulong x, long n, ulong p)
    2251             : {
    2252      296078 :   return Fl_powers_pre(x, n, p, (p & HIGHMASK)? get_Fl_red(p): 0);
    2253             : }
    2254             : 
    2255             : /**********************************************************************
    2256             :  **                    Powering  over (Z/NZ)^*, large N              **
    2257             :  **********************************************************************/
    2258             : typedef struct muldata {
    2259             :   GEN (*sqr)(void * E, GEN x);
    2260             :   GEN (*mul)(void * E, GEN x, GEN y);
    2261             :   GEN (*mul2)(void * E, GEN x);
    2262             : } muldata;
    2263             : 
    2264             : /* modified Barrett reduction with one fold */
    2265             : /* See Fast Modular Reduction, W. Hasenplaugh, G. Gaubatz, V. Gopal, ARITH 18 */
    2266             : 
    2267             : static GEN
    2268       15115 : Fp_invmBarrett(GEN p, long s)
    2269             : {
    2270       15115 :   GEN R, Q = dvmdii(int2n(3*s),p,&R);
    2271       15115 :   return mkvec2(Q,R);
    2272             : }
    2273             : 
    2274             : /* a <= (N-1)^2, 2^(2s-2) <= N < 2^(2s). Return 0 <= r < N such that
    2275             :  * a = r (mod N) */
    2276             : static GEN
    2277     9212010 : Fp_rem_mBarrett(GEN a, GEN B, long s, GEN N)
    2278             : {
    2279     9212010 :   pari_sp av = avma;
    2280     9212010 :   GEN P = gel(B, 1), Q = gel(B, 2); /* 2^(3s) = P N + Q, 0 <= Q < N */
    2281     9212010 :   long t = expi(P)+1; /* 2^(t-1) <= P < 2^t */
    2282     9212010 :   GEN u = shifti(a, -3*s), v = remi2n(a, 3*s); /* a = 2^(3s)u + v */
    2283     9212010 :   GEN A = addii(v, mulii(Q,u)); /* 0 <= A < 2^(3s+1) */
    2284     9212010 :   GEN q = shifti(mulii(shifti(A, t-3*s), P), -t); /* A/N - 4 < q <= A/N */
    2285     9212010 :   GEN r = subii(A, mulii(q, N));
    2286     9212010 :   GEN sr= subii(r,N);     /* 0 <= r < 4*N */
    2287     9212010 :   if (signe(sr)<0) return gerepileuptoint(av, r);
    2288     5027924 :   r=sr; sr = subii(r,N);  /* 0 <= r < 3*N */
    2289     5027924 :   if (signe(sr)<0) return gerepileuptoint(av, r);
    2290      186160 :   r=sr; sr = subii(r,N);  /* 0 <= r < 2*N */
    2291      186160 :   return gerepileuptoint(av, signe(sr)>=0 ? sr:r);
    2292             : }
    2293             : 
    2294             : /* Montgomery reduction */
    2295             : 
    2296             : INLINE ulong
    2297      670145 : init_montdata(GEN N) { return (ulong) -invmod2BIL(mod2BIL(N)); }
    2298             : 
    2299             : struct montred
    2300             : {
    2301             :   GEN N;
    2302             :   ulong inv;
    2303             : };
    2304             : 
    2305             : /* Montgomery reduction */
    2306             : static GEN
    2307    67879608 : _sqr_montred(void * E, GEN x)
    2308             : {
    2309    67879608 :   struct montred * D = (struct montred *) E;
    2310    67879608 :   return red_montgomery(sqri(x), D->N, D->inv);
    2311             : }
    2312             : 
    2313             : /* Montgomery reduction */
    2314             : static GEN
    2315     6958038 : _mul_montred(void * E, GEN x, GEN y)
    2316             : {
    2317     6958038 :   struct montred * D = (struct montred *) E;
    2318     6958038 :   return red_montgomery(mulii(x, y), D->N, D->inv);
    2319             : }
    2320             : 
    2321             : static GEN
    2322    11045744 : _mul2_montred(void * E, GEN x)
    2323             : {
    2324    11045744 :   struct montred * D = (struct montred *) E;
    2325    11045744 :   GEN z = shifti(_sqr_montred(E, x), 1);
    2326    11042728 :   long l = lgefint(D->N);
    2327    11682758 :   while (lgefint(z) > l) z = subii(z, D->N);
    2328    11043174 :   return z;
    2329             : }
    2330             : 
    2331             : static GEN
    2332    23030928 : _sqr_remii(void* N, GEN x)
    2333    23030928 : { return remii(sqri(x), (GEN) N); }
    2334             : 
    2335             : static GEN
    2336     1508958 : _mul_remii(void* N, GEN x, GEN y)
    2337     1508958 : { return remii(mulii(x, y), (GEN) N); }
    2338             : 
    2339             : static GEN
    2340     3175511 : _mul2_remii(void* N, GEN x)
    2341     3175511 : { return Fp_double(_sqr_remii(N, x), (GEN)N); }
    2342             : 
    2343             : struct redbarrett
    2344             : {
    2345             :   GEN iM, N;
    2346             :   long s;
    2347             : };
    2348             : 
    2349             : static GEN
    2350     8437640 : _sqr_remiibar(void *E, GEN x)
    2351             : {
    2352     8437640 :   struct redbarrett * D = (struct redbarrett *) E;
    2353     8437640 :   return Fp_rem_mBarrett(sqri(x), D->iM, D->s, D->N);
    2354             : }
    2355             : 
    2356             : static GEN
    2357      774370 : _mul_remiibar(void *E, GEN x, GEN y)
    2358             : {
    2359      774370 :   struct redbarrett * D = (struct redbarrett *) E;
    2360      774370 :   return Fp_rem_mBarrett(mulii(x, y), D->iM, D->s, D->N);
    2361             : }
    2362             : 
    2363             : static GEN
    2364     2079326 : _mul2_remiibar(void *E, GEN x)
    2365             : {
    2366     2079326 :   struct redbarrett * D = (struct redbarrett *) E;
    2367     2079326 :   return Fp_double(_sqr_remiibar(E, x), D->N);
    2368             : }
    2369             : 
    2370             : static long
    2371      864175 : Fp_select_red(GEN *y, ulong k, GEN N, long lN, muldata *D, void **pt_E)
    2372             : {
    2373      864175 :   if (lN >= Fp_POW_BARRETT_LIMIT && (k==0 || ((double)k)*expi(*y) > 2 + expi(N)))
    2374             :   {
    2375       15115 :     struct redbarrett * E = (struct redbarrett *) stack_malloc(sizeof(struct redbarrett));
    2376       15115 :     D->sqr = &_sqr_remiibar;
    2377       15115 :     D->mul = &_mul_remiibar;
    2378       15115 :     D->mul2 = &_mul2_remiibar;
    2379       15115 :     E->N = N;
    2380       15115 :     E->s = 1+(expi(N)>>1);
    2381       15115 :     E->iM = Fp_invmBarrett(N, E->s);
    2382       15115 :     *pt_E = (void*) E;
    2383       15115 :     return 0;
    2384             :   }
    2385      849060 :   else if (mod2(N) && lN < Fp_POW_REDC_LIMIT)
    2386             :   {
    2387      670140 :     struct montred * E = (struct montred *) stack_malloc(sizeof(struct montred));
    2388      670140 :     *y = remii(shifti(*y, bit_accuracy(lN)), N);
    2389      670145 :     D->sqr = &_sqr_montred;
    2390      670145 :     D->mul = &_mul_montred;
    2391      670145 :     D->mul2 = &_mul2_montred;
    2392      670145 :     E->N = N;
    2393      670145 :     E->inv = init_montdata(N);
    2394      670145 :     *pt_E = (void*) E;
    2395      670145 :     return 1;
    2396             :   }
    2397             :   else
    2398             :   {
    2399      178925 :     D->sqr = &_sqr_remii;
    2400      178925 :     D->mul = &_mul_remii;
    2401      178925 :     D->mul2 = &_mul2_remii;
    2402      178925 :     *pt_E = (void*) N;
    2403      178925 :     return 0;
    2404             :   }
    2405             : }
    2406             : 
    2407             : GEN
    2408     1752730 : Fp_powu(GEN A, ulong k, GEN N)
    2409             : {
    2410     1752730 :   long lN = lgefint(N);
    2411             :   int base_is_2, use_montgomery;
    2412             :   muldata D;
    2413             :   void *E;
    2414             :   pari_sp av;
    2415             : 
    2416     1752730 :   if (lN == 3) {
    2417      309653 :     ulong n = uel(N,2);
    2418      309653 :     return utoi( Fl_powu(umodiu(A, n), k, n) );
    2419             :   }
    2420     1443077 :   if (k <= 2)
    2421             :   { /* frequent special cases */
    2422      846617 :     if (k == 2) return Fp_sqr(A,N);
    2423      294912 :     if (k == 1) return A;
    2424           0 :     if (k == 0) return gen_1;
    2425             :   }
    2426      596460 :   av = avma; A = modii(A,N);
    2427      596460 :   base_is_2 = 0;
    2428      596460 :   if (lgefint(A) == 3) switch(A[2])
    2429             :   {
    2430        1055 :     case 1: set_avma(av); return gen_1;
    2431       34111 :     case 2:  base_is_2 = 1; break;
    2432             :   }
    2433             : 
    2434             :   /* TODO: Move this out of here and use for general modular computations */
    2435      595405 :   use_montgomery = Fp_select_red(&A, k, N, lN, &D, &E);
    2436      595405 :   if (base_is_2)
    2437       34111 :     A = gen_powu_fold_i(A, k, E, D.sqr, D.mul2);
    2438             :   else
    2439      561294 :     A = gen_powu_i(A, k, E, D.sqr, D.mul);
    2440      595405 :   if (use_montgomery)
    2441             :   {
    2442      499769 :     A = red_montgomery(A, N, ((struct montred *) E)->inv);
    2443      499769 :     if (cmpii(A, N) >= 0) A = subii(A,N);
    2444             :   }
    2445      595405 :   return gerepileuptoint(av, A);
    2446             : }
    2447             : 
    2448             : GEN
    2449     1346295 : Fp_pows(GEN A, long k, GEN N)
    2450             : {
    2451     1346295 :   if (lgefint(N) == 3) {
    2452     1322418 :     ulong n = N[2];
    2453     1322418 :     ulong a = umodiu(A, n);
    2454     1322419 :     if (k < 0) {
    2455       58634 :       a = Fl_inv(a, n);
    2456       58634 :       k = -k;
    2457             :     }
    2458     1322419 :     return utoi( Fl_powu(a, (ulong)k, n) );
    2459             :   }
    2460       23877 :   if (k < 0) { A = Fp_inv(A, N); k = -k; };
    2461       23877 :   return Fp_powu(A, (ulong)k, N);
    2462             : }
    2463             : 
    2464             : /* A^K mod N */
    2465             : GEN
    2466    36201863 : Fp_pow(GEN A, GEN K, GEN N)
    2467             : {
    2468             :   pari_sp av;
    2469    36201863 :   long s, lN = lgefint(N), sA, sy;
    2470             :   int base_is_2, use_montgomery;
    2471             :   GEN y;
    2472             :   muldata D;
    2473             :   void *E;
    2474             : 
    2475    36201863 :   s = signe(K);
    2476    36201863 :   if (!s) return dvdii(A,N)? gen_0: gen_1;
    2477    35171184 :   if (lN == 3 && lgefint(K) == 3)
    2478             :   {
    2479    34457800 :     ulong n = N[2], a = umodiu(A, n);
    2480    34458095 :     if (s < 0) a = Fl_inv(a, n);
    2481    34458222 :     if (a <= 1) return utoi(a); /* 0 or 1 */
    2482    30924489 :     return utoi(Fl_powu(a, uel(K,2), n));
    2483             :   }
    2484             : 
    2485      713384 :   av = avma;
    2486      713384 :   if (s < 0) y = Fp_inv(A,N);
    2487             :   else
    2488             :   {
    2489      711436 :     y = modii(A,N);
    2490      711625 :     if (!signe(y)) { set_avma(av); return gen_0; }
    2491             :   }
    2492      713573 :   if (lgefint(K) == 3) return gerepileuptoint(av, Fp_powu(y, K[2], N));
    2493             : 
    2494      268967 :   base_is_2 = 0;
    2495      268967 :   sy = abscmpii(y, shifti(N,-1)) > 0;
    2496      268968 :   if (sy) y = subii(N,y);
    2497      268971 :   sA = sy && mod2(K);
    2498      268971 :   if (lgefint(y) == 3) switch(y[2])
    2499             :   {
    2500         207 :     case 1:  set_avma(av); return sA ? subis(N,1): gen_1;
    2501      145078 :     case 2:  base_is_2 = 1; break;
    2502             :   }
    2503             : 
    2504             :   /* TODO: Move this out of here and use for general modular computations */
    2505      268764 :   use_montgomery = Fp_select_red(&y, 0UL, N, lN, &D, &E);
    2506      268777 :   if (base_is_2)
    2507      145089 :     y = gen_pow_fold_i(y, K, E, D.sqr, D.mul2);
    2508             :   else
    2509      123688 :     y = gen_pow_i(y, K, E, D.sqr, D.mul);
    2510      268787 :   if (use_montgomery)
    2511             :   {
    2512      170378 :     y = red_montgomery(y, N, ((struct montred *) E)->inv);
    2513      170381 :     if (cmpii(y,N) >= 0) y = subii(y,N);
    2514             :   }
    2515      268790 :   if (sA) y = subii(N, y);
    2516      268788 :   return gerepileuptoint(av,y);
    2517             : }
    2518             : 
    2519             : static GEN
    2520    14129421 : _Fp_mul(void *E, GEN x, GEN y) { return Fp_mul(x,y,(GEN)E); }
    2521             : static GEN
    2522     8134253 : _Fp_sqr(void *E, GEN x) { return Fp_sqr(x,(GEN)E); }
    2523             : static GEN
    2524       47162 : _Fp_one(void *E) { (void) E; return gen_1; }
    2525             : 
    2526             : GEN
    2527         105 : Fp_pow_init(GEN x, GEN n, long k, GEN p)
    2528         105 : { return gen_pow_init(x, n, k, (void*)p, &_Fp_sqr, &_Fp_mul); }
    2529             : 
    2530             : GEN
    2531       43694 : Fp_pow_table(GEN R, GEN n, GEN p)
    2532       43694 : { return gen_pow_table(R, n, (void*)p, &_Fp_one, &_Fp_mul); }
    2533             : 
    2534             : GEN
    2535        5931 : Fp_powers(GEN x, long n, GEN p)
    2536             : {
    2537        5931 :   if (lgefint(p) == 3)
    2538        2463 :     return Flv_to_ZV(Fl_powers(umodiu(x, uel(p, 2)), n, uel(p, 2)));
    2539        3468 :   return gen_powers(x, n, 1, (void*)p, _Fp_sqr, _Fp_mul, _Fp_one);
    2540             : }
    2541             : 
    2542             : GEN
    2543         497 : FpV_prod(GEN V, GEN p) { return gen_product(V, (void *)p, &_Fp_mul); }
    2544             : 
    2545             : static GEN
    2546    27869623 : _Fp_pow(void *E, GEN x, GEN n) { return Fp_pow(x,n,(GEN)E); }
    2547             : static GEN
    2548         153 : _Fp_rand(void *E) { return addiu(randomi(subiu((GEN)E,1)),1); }
    2549             : 
    2550             : static GEN Fp_easylog(void *E, GEN a, GEN g, GEN ord);
    2551             : static const struct bb_group Fp_star={_Fp_mul,_Fp_pow,_Fp_rand,hash_GEN,
    2552             :                                       equalii,equali1,Fp_easylog};
    2553             : 
    2554             : static GEN
    2555      889934 : _Fp_red(void *E, GEN x) { return Fp_red(x, (GEN)E); }
    2556             : static GEN
    2557     1175564 : _Fp_add(void *E, GEN x, GEN y) { (void) E; return addii(x,y); }
    2558             : static GEN
    2559     1086840 : _Fp_neg(void *E, GEN x) { (void) E; return negi(x); }
    2560             : static GEN
    2561      575348 : _Fp_rmul(void *E, GEN x, GEN y) { (void) E; return mulii(x,y); }
    2562             : static GEN
    2563       34307 : _Fp_inv(void *E, GEN x) { return Fp_inv(x,(GEN)E); }
    2564             : static int
    2565      260756 : _Fp_equal0(GEN x) { return signe(x)==0; }
    2566             : static GEN
    2567       19083 : _Fp_s(void *E, long x) { (void) E; return stoi(x); }
    2568             : 
    2569             : static const struct bb_field Fp_field={_Fp_red,_Fp_add,_Fp_rmul,_Fp_neg,
    2570             :                                         _Fp_inv,_Fp_equal0,_Fp_s};
    2571             : 
    2572        6963 : const struct bb_field *get_Fp_field(void **E, GEN p)
    2573        6963 : { *E = (void*)p; return &Fp_field; }
    2574             : 
    2575             : /*********************************************************************/
    2576             : /**               ORDER of INTEGERMOD x  in  (Z/nZ)*                **/
    2577             : /*********************************************************************/
    2578             : ulong
    2579      542566 : Fl_order(ulong a, ulong o, ulong p)
    2580             : {
    2581      542566 :   pari_sp av = avma;
    2582             :   GEN m, P, E;
    2583             :   long i;
    2584      542566 :   if (a==1) return 1;
    2585      445068 :   if (!o) o = p-1;
    2586      445068 :   m = factoru(o);
    2587      445068 :   P = gel(m,1);
    2588      445068 :   E = gel(m,2);
    2589     1265092 :   for (i = lg(P)-1; i; i--)
    2590             :   {
    2591      820024 :     ulong j, l = P[i], e = E[i], t = o / upowuu(l,e), y = Fl_powu(a, t, p);
    2592      820024 :     if (y == 1) o = t;
    2593      780406 :     else for (j = 1; j < e; j++)
    2594             :     {
    2595      386772 :       y = Fl_powu(y, l, p);
    2596      386772 :       if (y == 1) { o = t *  upowuu(l, j); break; }
    2597             :     }
    2598             :   }
    2599      445068 :   return gc_ulong(av, o);
    2600             : }
    2601             : 
    2602             : /*Find the exact order of a assuming a^o==1*/
    2603             : GEN
    2604      133460 : Fp_order(GEN a, GEN o, GEN p) {
    2605      133460 :   if (lgefint(p) == 3 && (!o || typ(o) == t_INT))
    2606             :   {
    2607       59273 :     ulong pp = p[2], oo = (o && lgefint(o)==3)? uel(o,2): pp-1;
    2608       59273 :     return utoi( Fl_order(umodiu(a, pp), oo, pp) );
    2609             :   }
    2610       74187 :   return gen_order(a, o, (void*)p, &Fp_star);
    2611             : }
    2612             : GEN
    2613          70 : Fp_factored_order(GEN a, GEN o, GEN p)
    2614          70 : { return gen_factored_order(a, o, (void*)p, &Fp_star); }
    2615             : 
    2616             : /* return order of a mod p^e, e > 0, pe = p^e */
    2617             : static GEN
    2618          70 : Zp_order(GEN a, GEN p, long e, GEN pe)
    2619             : {
    2620             :   GEN ap, op;
    2621          70 :   if (absequaliu(p, 2))
    2622             :   {
    2623          56 :     if (e == 1) return gen_1;
    2624          56 :     if (e == 2) return mod4(a) == 1? gen_1: gen_2;
    2625          49 :     if (mod4(a) == 1) op = gen_1; else { op = gen_2; a = Fp_sqr(a, pe); }
    2626             :   } else {
    2627          14 :     ap = (e == 1)? a: remii(a,p);
    2628          14 :     op = Fp_order(ap, subiu(p,1), p);
    2629          14 :     if (e == 1) return op;
    2630           0 :     a = Fp_pow(a, op, pe); /* 1 mod p */
    2631             :   }
    2632          49 :   if (equali1(a)) return op;
    2633           7 :   return mulii(op, powiu(p, e - Z_pval(subiu(a,1), p)));
    2634             : }
    2635             : 
    2636             : GEN
    2637          63 : znorder(GEN x, GEN o)
    2638             : {
    2639          63 :   pari_sp av = avma;
    2640             :   GEN b, a;
    2641             : 
    2642          63 :   if (typ(x) != t_INTMOD) pari_err_TYPE("znorder [t_INTMOD expected]",x);
    2643          56 :   b = gel(x,1); a = gel(x,2);
    2644          56 :   if (!equali1(gcdii(a,b))) pari_err_COPRIME("znorder", a,b);
    2645          49 :   if (!o)
    2646             :   {
    2647          35 :     GEN fa = Z_factor(b), P = gel(fa,1), E = gel(fa,2);
    2648          35 :     long i, l = lg(P);
    2649          35 :     o = gen_1;
    2650          70 :     for (i = 1; i < l; i++)
    2651             :     {
    2652          35 :       GEN p = gel(P,i);
    2653          35 :       long e = itos(gel(E,i));
    2654             : 
    2655          35 :       if (l == 2)
    2656          35 :         o = Zp_order(a, p, e, b);
    2657             :       else {
    2658           0 :         GEN pe = powiu(p,e);
    2659           0 :         o = lcmii(o, Zp_order(remii(a,pe), p, e, pe));
    2660             :       }
    2661             :     }
    2662          35 :     return gerepileuptoint(av, o);
    2663             :   }
    2664          14 :   return Fp_order(a, o, b);
    2665             : }
    2666             : 
    2667             : /*********************************************************************/
    2668             : /**               DISCRETE LOGARITHM  in  (Z/nZ)*                   **/
    2669             : /*********************************************************************/
    2670             : static GEN
    2671       56028 : Fp_log_halfgcd(ulong bnd, GEN C, GEN g, GEN p)
    2672             : {
    2673       56028 :   pari_sp av = avma;
    2674             :   GEN h1, h2, F, G;
    2675       56028 :   if (!Fp_ratlift(g,p,C,shifti(C,-1),&h1,&h2)) return gc_NULL(av);
    2676       33680 :   if ((F = Z_issmooth_fact(h1, bnd)) && (G = Z_issmooth_fact(h2, bnd)))
    2677             :   {
    2678         126 :     GEN M = cgetg(3, t_MAT);
    2679         126 :     gel(M,1) = vecsmall_concat(gel(F, 1),gel(G, 1));
    2680         126 :     gel(M,2) = vecsmall_concat(gel(F, 2),zv_neg_inplace(gel(G, 2)));
    2681         126 :     return gerepileupto(av, M);
    2682             :   }
    2683       33554 :   return gc_NULL(av);
    2684             : }
    2685             : 
    2686             : static GEN
    2687       56028 : Fp_log_find_rel(GEN b, ulong bnd, GEN C, GEN p, GEN *g, long *e)
    2688             : {
    2689             :   GEN rel;
    2690       56028 :   do { (*e)++; *g = Fp_mul(*g, b, p); rel = Fp_log_halfgcd(bnd, C, *g, p); }
    2691       56028 :   while (!rel);
    2692         126 :   return rel;
    2693             : }
    2694             : 
    2695             : struct Fp_log_rel
    2696             : {
    2697             :   GEN rel;
    2698             :   ulong prmax;
    2699             :   long nbrel, nbmax, nbgen;
    2700             : };
    2701             : 
    2702             : static long
    2703       59731 : tr(long i) { return odd(i) ? (i+1)>>1: -(i>>1); }
    2704             : 
    2705             : static long
    2706      169813 : rt(long i) { return i>0 ? 2*i-1: -2*i; }
    2707             : 
    2708             : /* add u^e */
    2709             : static void
    2710        2163 : addifsmooth1(struct Fp_log_rel *r, GEN z, long u, long e)
    2711             : {
    2712        2163 :   pari_sp av = avma;
    2713        2163 :   long off = r->prmax+1;
    2714        2163 :   GEN F = cgetg(3, t_MAT);
    2715        2163 :   gel(F,1) = vecsmall_append(gel(z,1), off+rt(u));
    2716        2163 :   gel(F,2) = vecsmall_append(gel(z,2), e);
    2717        2163 :   gel(r->rel,++r->nbrel) = gerepileupto(av, F);
    2718        2163 : }
    2719             : 
    2720             : /* add u^-1 v^-1 */
    2721             : static void
    2722       83825 : addifsmooth2(struct Fp_log_rel *r, GEN z, long u, long v)
    2723             : {
    2724       83825 :   pari_sp av = avma;
    2725       83825 :   long off = r->prmax+1;
    2726       83825 :   GEN P = mkvecsmall2(off+rt(u),off+rt(v)), E = mkvecsmall2(-1,-1);
    2727       83825 :   GEN F = cgetg(3, t_MAT);
    2728       83825 :   gel(F,1) = vecsmall_concat(gel(z,1), P);
    2729       83825 :   gel(F,2) = vecsmall_concat(gel(z,2), E);
    2730       83825 :   gel(r->rel,++r->nbrel) = gerepileupto(av, F);
    2731       83825 : }
    2732             : 
    2733             : /* Let p=C^2+c
    2734             :  * Solve h = (C+x)*(C+a)-p = 0 [mod l]
    2735             :  * h= -c+x*(C+a)+C*a = 0  [mod l]
    2736             :  * x = (c-C*a)/(C+a) [mod l]
    2737             :  * h = -c+C*(x+a)+a*x */
    2738             : GEN
    2739       30249 : Fp_log_sieve_worker(long a, long prmax, GEN C, GEN c, GEN Ci, GEN ci, GEN pi, GEN sz)
    2740             : {
    2741       30249 :   pari_sp ltop = avma;
    2742       30249 :   long i, j, th, n = lg(pi)-1, rel = 1, ab = labs(a), ae;
    2743       30249 :   GEN sieve = zero_zv(2*ab+2)+1+ab;
    2744       30259 :   GEN L = cgetg(1+2*ab+2, t_VEC);
    2745       30252 :   pari_sp av = avma;
    2746       30252 :   GEN z, h = addis(C,a);
    2747       30253 :   if ((z = Z_issmooth_fact(h, prmax)))
    2748             :   {
    2749        2169 :     gel(L, rel++) = mkvec2(z, mkvecsmall3(1, a, -1));
    2750        2169 :     av = avma;
    2751             :   }
    2752    12475925 :   for (i=1; i<=n; i++)
    2753             :   {
    2754    12447544 :     ulong li = pi[i], s = sz[i], al = smodss(a,li);
    2755    12435795 :     ulong iv = Fl_invsafe(Fl_add(Ci[i],al,li),li);
    2756             :     long u;
    2757    12695548 :     if (!iv) continue;
    2758    12381020 :     u = Fl_add(Fl_mul(Fl_sub(ci[i],Fl_mul(Ci[i],al,li),li), iv ,li),ab%li,li)-ab;
    2759    46256952 :     for(j = u; j<=ab; j+=li) sieve[j] += s;
    2760             :   }
    2761       28381 :   if (a)
    2762             :   {
    2763       30136 :     long e = expi(mulis(C,a));
    2764       30161 :     th = e - expu(e) - 1;
    2765          54 :   } else th = -1;
    2766       30251 :   ae = a>=0 ? ab-1: ab;
    2767    15516109 :   for (j = 1-ab; j <= ae; j++)
    2768    15484794 :     if (sieve[j]>=th)
    2769             :     {
    2770      108859 :       GEN h = absi(addis(subii(mulis(C,a+j),c), a*j));
    2771      108713 :       if ((z = Z_issmooth_fact(h, prmax)))
    2772             :       {
    2773      106473 :         gel(L, rel++) = mkvec2(z, mkvecsmall3(2, a, j));
    2774      106502 :         av = avma;
    2775        2293 :       } else set_avma(av);
    2776             :     }
    2777             :   /* j = a */
    2778       31315 :   if (sieve[a]>=th)
    2779             :   {
    2780         448 :     GEN h = absi(addiu(subii(mulis(C,2*a),c), a*a));
    2781         448 :     if ((z = Z_issmooth_fact(h, prmax)))
    2782         364 :       gel(L, rel++) = mkvec2(z, mkvecsmall3(1, a, -2));
    2783             :   }
    2784       31315 :   setlg(L, rel); return gerepilecopy(ltop, L);
    2785             : }
    2786             : 
    2787             : static long
    2788          63 : Fp_log_sieve(struct Fp_log_rel *r, GEN C, GEN c, GEN Ci, GEN ci, GEN pi, GEN sz)
    2789             : {
    2790             :   struct pari_mt pt;
    2791          63 :   long i, j, nb = 0;
    2792          63 :   GEN worker = snm_closure(is_entry("_Fp_log_sieve_worker"),
    2793             :                mkvecn(7, utoi(r->prmax), C, c, Ci, ci, pi, sz));
    2794          63 :   long running, pending = 0;
    2795          63 :   GEN W = zerovec(r->nbgen);
    2796          63 :   mt_queue_start_lim(&pt, worker, r->nbgen);
    2797       30459 :   for (i = 0; (running = (i < r->nbgen)) || pending; i++)
    2798             :   {
    2799             :     GEN done;
    2800             :     long idx;
    2801       30396 :     mt_queue_submit(&pt, i, running ? mkvec(stoi(tr(i))): NULL);
    2802       30396 :     done = mt_queue_get(&pt, &idx, &pending);
    2803       30396 :     if (!done || lg(done)==1) continue;
    2804       27636 :     gel(W, idx+1) = done;
    2805       27636 :     nb += lg(done)-1;
    2806       27636 :     if (DEBUGLEVEL && (i&127)==0)
    2807           0 :       err_printf("%ld%% ",100*nb/r->nbmax);
    2808             :   }
    2809          63 :   mt_queue_end(&pt);
    2810       26362 :   for(j = 1; j <= r->nbgen && r->nbrel < r->nbmax; j++)
    2811             :   {
    2812             :     long ll, m;
    2813       26299 :     GEN L = gel(W,j);
    2814       26299 :     if (isintzero(L)) continue;
    2815       23681 :     ll = lg(L);
    2816      109669 :     for (m=1; m<ll && r->nbrel < r->nbmax ; m++)
    2817             :     {
    2818       85988 :       GEN Lm = gel(L,m), h = gel(Lm, 1), v = gel(Lm, 2);
    2819       85988 :       if (v[1] == 1)
    2820        2163 :         addifsmooth1(r, h, v[2], v[3]);
    2821             :       else
    2822       83825 :         addifsmooth2(r, h, v[2], v[3]);
    2823             :     }
    2824             :   }
    2825          63 :   return j;
    2826             : }
    2827             : 
    2828             : static GEN
    2829         837 : ECP_psi(GEN x, GEN y)
    2830             : {
    2831         837 :   long prec = realprec(x);
    2832         837 :   GEN lx = glog(x, prec), ly = glog(y, prec);
    2833         837 :   GEN u = gdiv(lx, ly);
    2834         837 :   return gpow(u, gneg(u),prec);
    2835             : }
    2836             : 
    2837             : struct computeG
    2838             : {
    2839             :   GEN C;
    2840             :   long bnd, nbi;
    2841             : };
    2842             : 
    2843             : static GEN
    2844         837 : _computeG(void *E, GEN gen)
    2845             : {
    2846         837 :   struct computeG * d = (struct computeG *) E;
    2847         837 :   GEN ps = ECP_psi(gmul(gen,d->C), stoi(d->bnd));
    2848         837 :   return gsub(gmul(gsqr(gen),ps),gmulgs(gaddgs(gen,d->nbi),3));
    2849             : }
    2850             : 
    2851             : static long
    2852          63 : compute_nbgen(GEN C, long bnd, long nbi)
    2853             : {
    2854             :   struct computeG d;
    2855          63 :   d.C = shifti(C, 1);
    2856          63 :   d.bnd = bnd;
    2857          63 :   d.nbi = nbi;
    2858          63 :   return itos(ground(zbrent((void*)&d, _computeG, gen_2, stoi(bnd), DEFAULTPREC)));
    2859             : }
    2860             : 
    2861             : static GEN
    2862        1714 : _psi(void*E, GEN y)
    2863             : {
    2864        1714 :   GEN lx = (GEN) E;
    2865        1714 :   long prec = realprec(lx);
    2866        1714 :   GEN ly = glog(y, prec);
    2867        1714 :   GEN u = gdiv(lx, ly);
    2868        1714 :   return gsub(gdiv(y ,ly), gpow(u, u, prec));
    2869             : }
    2870             : 
    2871             : static GEN
    2872          63 : opt_param(GEN x, long prec)
    2873             : {
    2874          63 :   return zbrent((void*)glog(x,prec), _psi, gen_2, x, prec);
    2875             : }
    2876             : 
    2877             : static GEN
    2878          63 : check_kernel(long nbg, long N, long prmax, GEN C, GEN M, GEN p, GEN m)
    2879             : {
    2880          63 :   pari_sp av = avma;
    2881          63 :   long lM = lg(M)-1, nbcol = lM;
    2882          63 :   long tbs = maxss(1, expu(nbg/expi(m)));
    2883             :   for (;;)
    2884          42 :   {
    2885         105 :     GEN K = FpMs_leftkernel_elt_col(M, nbcol, N, m);
    2886             :     GEN tab;
    2887         105 :     long i, f=0;
    2888         105 :     long l = lg(K), lm = lgefint(m);
    2889         105 :     GEN idx = diviiexact(subiu(p,1),m), g;
    2890             :     pari_timer ti;
    2891         105 :     if (DEBUGLEVEL) timer_start(&ti);
    2892         210 :     for(i=1; i<l; i++)
    2893         210 :       if (signe(gel(K,i)))
    2894         105 :         break;
    2895         105 :     g = Fp_pow(utoi(i), idx, p);
    2896         105 :     tab = Fp_pow_init(g, p, tbs, p);
    2897         105 :     K = FpC_Fp_mul(K, Fp_inv(gel(K,i), m), m);
    2898      121520 :     for(i=1; i<l; i++)
    2899             :     {
    2900      121415 :       GEN k = gel(K,i);
    2901      121415 :       GEN j = i<=prmax ? utoi(i): addis(C,tr(i-(prmax+1)));
    2902      121415 :       if (signe(k)==0 || !equalii(Fp_pow_table(tab, k, p), Fp_pow(j, idx, p)))
    2903       82369 :         gel(K,i) = cgetineg(lm);
    2904             :       else
    2905       39046 :         f++;
    2906             :     }
    2907         105 :     if (DEBUGLEVEL) timer_printf(&ti,"found %ld/%ld logs", f, nbg);
    2908         105 :     if(f > (nbg>>1)) return gerepileupto(av, K);
    2909       10024 :     for(i=1; i<=nbcol; i++)
    2910             :     {
    2911        9982 :       long a = 1+random_Fl(lM);
    2912        9982 :       swap(gel(M,a),gel(M,i));
    2913             :     }
    2914          42 :     if (4*nbcol>5*nbg) nbcol = nbcol*9/10;
    2915             :   }
    2916             : }
    2917             : 
    2918             : static GEN
    2919         126 : Fp_log_find_ind(GEN a, GEN K, long prmax, GEN C, GEN p, GEN m)
    2920             : {
    2921         126 :   pari_sp av=avma;
    2922         126 :   GEN aa = gen_1;
    2923         126 :   long AV = 0;
    2924             :   for(;;)
    2925           0 :   {
    2926         126 :     GEN A = Fp_log_find_rel(a, prmax, C, p, &aa, &AV);
    2927         126 :     GEN F = gel(A,1), E = gel(A,2);
    2928         126 :     GEN Ao = gen_0;
    2929         126 :     long i, l = lg(F);
    2930         807 :     for(i=1; i<l; i++)
    2931             :     {
    2932         681 :       GEN Ki = gel(K,F[i]);
    2933         681 :       if (signe(Ki)<0) break;
    2934         681 :       Ao = addii(Ao, mulis(Ki, E[i]));
    2935             :     }
    2936         126 :     if (i==l) return Fp_divu(Ao, AV, m);
    2937           0 :     aa = gerepileuptoint(av, aa);
    2938             :   }
    2939             : }
    2940             : 
    2941             : static GEN
    2942          63 : Fp_log_index(GEN a, GEN b, GEN m, GEN p)
    2943             : {
    2944          63 :   pari_sp av = avma, av2;
    2945          63 :   long i, j, nbi, nbr = 0, nbrow, nbg;
    2946             :   GEN C, c, Ci, ci, pi, pr, sz, l, Ao, Bo, K, d, p_1;
    2947             :   pari_timer ti;
    2948             :   struct Fp_log_rel r;
    2949          63 :   ulong bnds = itou(roundr_safe(opt_param(sqrti(p),DEFAULTPREC)));
    2950          63 :   ulong bnd = 4*bnds;
    2951          63 :   if (!bnds || cmpii(sqru(bnds),m)>=0) return NULL;
    2952             : 
    2953          63 :   p_1 = subiu(p,1);
    2954          63 :   if (!is_pm1(gcdii(m,diviiexact(p_1,m))))
    2955           0 :     m = diviiexact(p_1, Z_ppo(p_1, m));
    2956          63 :   pr = primes_upto_zv(bnd);
    2957          63 :   nbi = lg(pr)-1;
    2958          63 :   C = sqrtremi(p, &c);
    2959          63 :   av2 = avma;
    2960       12796 :   for (i = 1; i <= nbi; ++i)
    2961             :   {
    2962       12733 :     ulong lp = pr[i];
    2963       26894 :     while (lp <= bnd)
    2964             :     {
    2965       14161 :       nbr++;
    2966       14161 :       lp *= pr[i];
    2967             :     }
    2968             :   }
    2969          63 :   pi = cgetg(nbr+1,t_VECSMALL);
    2970          63 :   Ci = cgetg(nbr+1,t_VECSMALL);
    2971          63 :   ci = cgetg(nbr+1,t_VECSMALL);
    2972          63 :   sz = cgetg(nbr+1,t_VECSMALL);
    2973       12796 :   for (i = 1, j = 1; i <= nbi; ++i)
    2974             :   {
    2975       12733 :     ulong lp = pr[i], sp = expu(2*lp-1);
    2976       26894 :     while (lp <= bnd)
    2977             :     {
    2978       14161 :       pi[j] = lp;
    2979       14161 :       Ci[j] = umodiu(C, lp);
    2980       14161 :       ci[j] = umodiu(c, lp);
    2981       14161 :       sz[j] = sp;
    2982       14161 :       lp *= pr[i];
    2983       14161 :       j++;
    2984             :     }
    2985             :   }
    2986          63 :   r.nbrel = 0;
    2987          63 :   r.nbgen = compute_nbgen(C, bnd, nbi);
    2988          63 :   r.nbmax = 2*(nbi+r.nbgen);
    2989          63 :   r.rel = cgetg(r.nbmax+1,t_VEC);
    2990          63 :   r.prmax = pr[nbi];
    2991          63 :   if (DEBUGLEVEL)
    2992             :   {
    2993           0 :     err_printf("bnd=%lu Size FB=%ld extra gen=%ld \n", bnd, nbi, r.nbgen);
    2994           0 :     timer_start(&ti);
    2995             :   }
    2996          63 :   nbg = Fp_log_sieve(&r, C, c, Ci, ci, pi, sz);
    2997          63 :   nbrow = r.prmax + nbg;
    2998          63 :   if (DEBUGLEVEL)
    2999             :   {
    3000           0 :     err_printf("\n");
    3001           0 :     timer_printf(&ti," %ld relations, %ld generators", r.nbrel, nbi+nbg);
    3002             :   }
    3003          63 :   setlg(r.rel,r.nbrel+1);
    3004          63 :   r.rel = gerepilecopy(av2, r.rel);
    3005          63 :   K = check_kernel(nbi+nbrow-r.prmax, nbrow, r.prmax, C, r.rel, p, m);
    3006          63 :   if (DEBUGLEVEL) timer_start(&ti);
    3007          63 :   Ao = Fp_log_find_ind(a, K, r.prmax, C, p, m);
    3008          63 :   if (DEBUGLEVEL) timer_printf(&ti," log element");
    3009          63 :   Bo = Fp_log_find_ind(b, K, r.prmax, C, p, m);
    3010          63 :   if (DEBUGLEVEL) timer_printf(&ti," log generator");
    3011          63 :   d = gcdii(Ao,Bo);
    3012          63 :   l = Fp_div(diviiexact(Ao, d) ,diviiexact(Bo, d), m);
    3013          63 :   if (!equalii(a,Fp_pow(b,l,p))) pari_err_BUG("Fp_log_index");
    3014          63 :   return gerepileuptoint(av, l);
    3015             : }
    3016             : 
    3017             : static int
    3018     4664059 : Fp_log_use_index(long e, long p)
    3019             : {
    3020     4664059 :   return (e >= 27 && 20*(p+6)<=e*e);
    3021             : }
    3022             : 
    3023             : /* Trivial cases a = 1, -1. Return x s.t. g^x = a or [] if no such x exist */
    3024             : static GEN
    3025     8464017 : Fp_easylog(void *E, GEN a, GEN g, GEN ord)
    3026             : {
    3027     8464017 :   pari_sp av = avma;
    3028     8464017 :   GEN p = (GEN)E;
    3029             :   /* assume a reduced mod p, p not necessarily prime */
    3030     8464017 :   if (equali1(a)) return gen_0;
    3031             :   /* p > 2 */
    3032     5440672 :   if (equalii(subiu(p,1), a))  /* -1 */
    3033             :   {
    3034             :     pari_sp av2;
    3035             :     GEN t;
    3036     1323450 :     ord = get_arith_Z(ord);
    3037     1323450 :     if (mpodd(ord)) { set_avma(av); return cgetg(1, t_VEC); } /* no solution */
    3038     1323436 :     t = shifti(ord,-1); /* only possible solution */
    3039     1323436 :     av2 = avma;
    3040     1323436 :     if (!equalii(Fp_pow(g, t, p), a)) { set_avma(av); return cgetg(1, t_VEC); }
    3041     1323408 :     set_avma(av2); return gerepileuptoint(av, t);
    3042             :   }
    3043     4117227 :   if (typ(ord)==t_INT && BPSW_psp(p) && Fp_log_use_index(expi(ord),expi(p)))
    3044          63 :     return Fp_log_index(a, g, ord, p);
    3045     4117164 :   return gc_NULL(av); /* not easy */
    3046             : }
    3047             : 
    3048             : GEN
    3049     3926445 : Fp_log(GEN a, GEN g, GEN ord, GEN p)
    3050             : {
    3051     3926445 :   GEN v = get_arith_ZZM(ord);
    3052     3926412 :   GEN F = gmael(v,2,1);
    3053     3926412 :   long lF = lg(F)-1, lmax;
    3054     3926412 :   if (lF == 0) return equali1(a)? gen_0: cgetg(1, t_VEC);
    3055     3926384 :   lmax = expi(gel(F,lF));
    3056     3926381 :   if (BPSW_psp(p) && Fp_log_use_index(lmax,expi(p)))
    3057          91 :     v = mkvec2(gel(v,1),ZM_famat_limit(gel(v,2),int2n(27)));
    3058     3926374 :   return gen_PH_log(a,g,v,(void*)p,&Fp_star);
    3059             : }
    3060             : 
    3061             : /* assume !(p & HIGHMASK) */
    3062             : static ulong
    3063      132738 : Fl_log_naive(ulong a, ulong g, ulong ord, ulong p)
    3064             : {
    3065      132738 :   ulong i, h=1;
    3066      365021 :   for (i = 0; i < ord; i++, h = (h * g) % p)
    3067      365021 :     if (a==h) return i;
    3068           0 :   return ~0UL;
    3069             : }
    3070             : 
    3071             : static ulong
    3072       25146 : Fl_log_naive_pre(ulong a, ulong g, ulong ord, ulong p, ulong pi)
    3073             : {
    3074       25146 :   ulong i, h=1;
    3075       64430 :   for (i = 0; i < ord; i++, h = Fl_mul_pre(h, g, p, pi))
    3076       64430 :     if (a==h) return i;
    3077           0 :   return ~0UL;
    3078             : }
    3079             : 
    3080             : static ulong
    3081           0 : Fl_log_Fp(ulong a, ulong g, ulong ord, ulong p)
    3082             : {
    3083           0 :   pari_sp av = avma;
    3084           0 :   GEN r = Fp_log(utoi(a),utoi(g),utoi(ord),utoi(p));
    3085           0 :   return gc_ulong(av, typ(r)==t_INT ? itou(r): ~0UL);
    3086             : }
    3087             : 
    3088             : /* allow pi = 0 */
    3089             : ulong
    3090       25557 : Fl_log_pre(ulong a, ulong g, ulong ord, ulong p, ulong pi)
    3091             : {
    3092       25557 :   if (!pi) return Fl_log(a, g, ord, p);
    3093       25146 :   if (ord <= 200) return Fl_log_naive_pre(a, g, ord, p, pi);
    3094           0 :   return Fl_log_Fp(a, g, ord, p);
    3095             : }
    3096             : 
    3097             : ulong
    3098      132738 : Fl_log(ulong a, ulong g, ulong ord, ulong p)
    3099             : {
    3100      132738 :   if (ord <= 200)
    3101           0 :     return (p&HIGHMASK)? Fl_log_naive_pre(a, g, ord, p, get_Fl_red(p))
    3102      132738 :                        : Fl_log_naive(a, g, ord, p);
    3103           0 :   return Fl_log_Fp(a, g, ord, p);
    3104             : }
    3105             : 
    3106             : /* find x such that h = g^x mod N > 1, N = prod_{i <= l} P[i]^E[i], P[i] prime.
    3107             :  * PHI[l] = eulerphi(N / P[l]^E[l]).   Destroys P/E */
    3108             : static GEN
    3109         126 : znlog_rec(GEN h, GEN g, GEN N, GEN P, GEN E, GEN PHI)
    3110             : {
    3111         126 :   long l = lg(P) - 1, e = E[l];
    3112         126 :   GEN p = gel(P, l), phi = gel(PHI,l), pe = e == 1? p: powiu(p, e);
    3113             :   GEN a,b, hp,gp, hpe,gpe, ogpe; /* = order(g mod p^e) | p^(e-1)(p-1) */
    3114             : 
    3115         126 :   if (l == 1) {
    3116          98 :     hpe = h;
    3117          98 :     gpe = g;
    3118             :   } else {
    3119          28 :     hpe = modii(h, pe);
    3120          28 :     gpe = modii(g, pe);
    3121             :   }
    3122         126 :   if (e == 1) {
    3123          42 :     hp = hpe;
    3124          42 :     gp = gpe;
    3125             :   } else {
    3126          84 :     hp = remii(hpe, p);
    3127          84 :     gp = remii(gpe, p);
    3128             :   }
    3129         126 :   if (hp == gen_0 || gp == gen_0) return NULL;
    3130         105 :   if (absequaliu(p, 2))
    3131             :   {
    3132          35 :     GEN n = int2n(e);
    3133          35 :     ogpe = Zp_order(gpe, gen_2, e, n);
    3134          35 :     a = Fp_log(hpe, gpe, ogpe, n);
    3135          35 :     if (typ(a) != t_INT) return NULL;
    3136             :   }
    3137             :   else
    3138             :   { /* Avoid black box groups: (Z/p^2)^* / (Z/p)^* ~ (Z/pZ, +), where DL
    3139             :        is trivial */
    3140             :     /* [order(gp), factor(order(gp))] */
    3141          70 :     GEN v = Fp_factored_order(gp, subiu(p,1), p);
    3142          70 :     GEN ogp = gel(v,1);
    3143          70 :     if (!equali1(Fp_pow(hp, ogp, p))) return NULL;
    3144          70 :     a = Fp_log(hp, gp, v, p);
    3145          70 :     if (typ(a) != t_INT) return NULL;
    3146          70 :     if (e == 1) ogpe = ogp;
    3147             :     else
    3148             :     { /* find a s.t. g^a = h (mod p^e), p odd prime, e > 0, (h,p) = 1 */
    3149             :       /* use p-adic log: O(log p + e) mul*/
    3150             :       long vpogpe, vpohpe;
    3151             : 
    3152          28 :       hpe = Fp_mul(hpe, Fp_pow(gpe, negi(a), pe), pe);
    3153          28 :       gpe = Fp_pow(gpe, ogp, pe);
    3154             :       /* g,h = 1 mod p; compute b s.t. h = g^b */
    3155             : 
    3156             :       /* v_p(order g mod pe) */
    3157          28 :       vpogpe = equali1(gpe)? 0: e - Z_pval(subiu(gpe,1), p);
    3158             :       /* v_p(order h mod pe) */
    3159          28 :       vpohpe = equali1(hpe)? 0: e - Z_pval(subiu(hpe,1), p);
    3160          28 :       if (vpohpe > vpogpe) return NULL;
    3161             : 
    3162          28 :       ogpe = mulii(ogp, powiu(p, vpogpe)); /* order g mod p^e */
    3163          28 :       if (is_pm1(gpe)) return is_pm1(hpe)? a: NULL;
    3164          28 :       b = gdiv(Qp_log(cvtop(hpe, p, e)), Qp_log(cvtop(gpe, p, e)));
    3165          28 :       a = addii(a, mulii(ogp, padic_to_Q(b)));
    3166             :     }
    3167             :   }
    3168             :   /* gp^a = hp => x = a mod ogpe => generalized Pohlig-Hellman strategy */
    3169          91 :   if (l == 1) return a;
    3170             : 
    3171          28 :   N = diviiexact(N, pe); /* make N coprime to p */
    3172          28 :   h = Fp_mul(h, Fp_pow(g, modii(negi(a), phi), N), N);
    3173          28 :   g = Fp_pow(g, modii(ogpe, phi), N);
    3174          28 :   setlg(P, l); /* remove last element */
    3175          28 :   setlg(E, l);
    3176          28 :   b = znlog_rec(h, g, N, P, E, PHI);
    3177          28 :   if (!b) return NULL;
    3178          28 :   return addmulii(a, b, ogpe);
    3179             : }
    3180             : 
    3181             : static GEN
    3182          98 : get_PHI(GEN P, GEN E)
    3183             : {
    3184          98 :   long i, l = lg(P);
    3185          98 :   GEN PHI = cgetg(l, t_VEC);
    3186          98 :   gel(PHI,1) = gen_1;
    3187         126 :   for (i=1; i<l-1; i++)
    3188             :   {
    3189          28 :     GEN t, p = gel(P,i);
    3190          28 :     long e = E[i];
    3191          28 :     t = mulii(powiu(p, e-1), subiu(p,1));
    3192          28 :     if (i > 1) t = mulii(t, gel(PHI,i));
    3193          28 :     gel(PHI,i+1) = t;
    3194             :   }
    3195          98 :   return PHI;
    3196             : }
    3197             : 
    3198             : GEN
    3199         238 : znlog(GEN h, GEN g, GEN o)
    3200             : {
    3201         238 :   pari_sp av = avma;
    3202             :   GEN N, fa, P, E, x;
    3203         238 :   switch (typ(g))
    3204             :   {
    3205          28 :     case t_PADIC:
    3206             :     {
    3207          28 :       GEN p = padic_p(g);
    3208          28 :       long v = valp(g);
    3209          28 :       if (v < 0) pari_err_DIM("znlog");
    3210          28 :       if (v > 0) {
    3211           0 :         long k = gvaluation(h, p);
    3212           0 :         if (k % v) return cgetg(1,t_VEC);
    3213           0 :         k /= v;
    3214           0 :         if (!gequal(h, gpowgs(g,k))) { set_avma(av); return cgetg(1,t_VEC); }
    3215           0 :         return gc_stoi(av, k);
    3216             :       }
    3217          28 :       N = padic_pd(g);
    3218          28 :       g = Rg_to_Fp(g, N);
    3219          28 :       break;
    3220             :     }
    3221         203 :     case t_INTMOD:
    3222         203 :       N = gel(g,1);
    3223         203 :       g = gel(g,2); break;
    3224           7 :     default: pari_err_TYPE("znlog", g);
    3225             :       return NULL; /* LCOV_EXCL_LINE */
    3226             :   }
    3227         231 :   if (equali1(N)) { set_avma(av); return gen_0; }
    3228         231 :   h = Rg_to_Fp(h, N);
    3229         224 :   if (o) return gerepileupto(av, Fp_log(h, g, o, N));
    3230          98 :   fa = Z_factor(N);
    3231          98 :   P = gel(fa,1);
    3232          98 :   E = vec_to_vecsmall(gel(fa,2));
    3233          98 :   x = znlog_rec(h, g, N, P, E, get_PHI(P,E));
    3234          98 :   if (!x) { set_avma(av); return cgetg(1,t_VEC); }
    3235          63 :   return gerepileuptoint(av, x);
    3236             : }
    3237             : 
    3238             : GEN
    3239      173541 : Fp_sqrtn(GEN a, GEN n, GEN p, GEN *zeta)
    3240             : {
    3241      173541 :   if (lgefint(p)==3)
    3242             :   {
    3243      172921 :     long nn = itos_or_0(n);
    3244      172921 :     if (nn)
    3245             :     {
    3246      172921 :       ulong pp = p[2];
    3247             :       ulong uz;
    3248      172921 :       ulong r = Fl_sqrtn(umodiu(a,pp),nn,pp, zeta ? &uz:NULL);
    3249      172900 :       if (r==ULONG_MAX) return NULL;
    3250      172858 :       if (zeta) *zeta = utoi(uz);
    3251      172858 :       return utoi(r);
    3252             :     }
    3253             :   }
    3254         620 :   a = modii(a,p);
    3255         620 :   if (!signe(a))
    3256             :   {
    3257           0 :     if (zeta) *zeta = gen_1;
    3258           0 :     if (signe(n) < 0) pari_err_INV("Fp_sqrtn", mkintmod(gen_0,p));
    3259           0 :     return gen_0;
    3260             :   }
    3261         620 :   if (absequaliu(n,2))
    3262             :   {
    3263         420 :     if (zeta) *zeta = subiu(p,1);
    3264         420 :     return signe(n) > 0 ? Fp_sqrt(a,p): Fp_sqrt(Fp_inv(a, p),p);
    3265             :   }
    3266         200 :   return gen_Shanks_sqrtn(a,n,subiu(p,1),zeta,(void*)p,&Fp_star);
    3267             : }
    3268             : 
    3269             : /*********************************************************************/
    3270             : /**                              FACTORIAL                          **/
    3271             : /*********************************************************************/
    3272             : GEN
    3273       90640 : mulu_interval_step(ulong a, ulong b, ulong step)
    3274             : {
    3275       90640 :   pari_sp av = avma;
    3276             :   ulong k, l, N, n;
    3277             :   long lx;
    3278             :   GEN x;
    3279             : 
    3280       90640 :   if (!a) return gen_0;
    3281       90640 :   if (step == 1) return mulu_interval(a, b);
    3282       90640 :   n = 1 + (b-a) / step;
    3283       90640 :   b -= (b-a) % step;
    3284       90640 :   if (n < 61)
    3285             :   {
    3286       89256 :     if (n == 1) return utoipos(a);
    3287       68713 :     x = muluu(a,a+step); if (n == 2) return x;
    3288      539092 :     for (k=a+2*step; k<=b; k+=step) x = mului(k,x);
    3289       53922 :     return gerepileuptoint(av, x);
    3290             :   }
    3291             :   /* step | b-a */
    3292        1384 :   lx = 1; x = cgetg(2 + n/2, t_VEC);
    3293        1384 :   N = b + a;
    3294        1384 :   for (k = a;; k += step)
    3295             :   {
    3296      227455 :     l = N - k; if (l <= k) break;
    3297      226071 :     gel(x,lx++) = muluu(k,l);
    3298             :   }
    3299        1384 :   if (l == k) gel(x,lx++) = utoipos(k);
    3300        1384 :   setlg(x, lx);
    3301        1384 :   return gerepileuptoint(av, ZV_prod(x));
    3302             : }
    3303             : /* return a * (a+1) * ... * b. Assume a <= b  [ note: factoring out powers of 2
    3304             :  * first is slower ... ] */
    3305             : GEN
    3306      158926 : mulu_interval(ulong a, ulong b)
    3307             : {
    3308      158926 :   pari_sp av = avma;
    3309             :   ulong k, l, N, n;
    3310             :   long lx;
    3311             :   GEN x;
    3312             : 
    3313      158926 :   if (!a) return gen_0;
    3314      158926 :   n = b - a + 1;
    3315      158926 :   if (n < 61)
    3316             :   {
    3317      158196 :     if (n == 1) return utoipos(a);
    3318      107880 :     x = muluu(a,a+1); if (n == 2) return x;
    3319      403545 :     for (k=a+2; k<b; k++) x = mului(k,x);
    3320             :     /* avoid k <= b: broken if b = ULONG_MAX */
    3321       93787 :     return gerepileuptoint(av, mului(b,x));
    3322             :   }
    3323         730 :   lx = 1; x = cgetg(2 + n/2, t_VEC);
    3324         732 :   N = b + a;
    3325         732 :   for (k = a;; k++)
    3326             :   {
    3327       27595 :     l = N - k; if (l <= k) break;
    3328       26864 :     gel(x,lx++) = muluu(k,l);
    3329             :   }
    3330         731 :   if (l == k) gel(x,lx++) = utoipos(k);
    3331         731 :   setlg(x, lx);
    3332         731 :   return gerepileuptoint(av, ZV_prod(x));
    3333             : }
    3334             : GEN
    3335         560 : muls_interval(long a, long b)
    3336             : {
    3337         560 :   pari_sp av = avma;
    3338         560 :   long lx, k, l, N, n = b - a + 1;
    3339             :   GEN x;
    3340             : 
    3341         560 :   if (a <= 0 && b >= 0) return gen_0;
    3342         287 :   if (n < 61)
    3343             :   {
    3344         287 :     x = stoi(a);
    3345         511 :     for (k=a+1; k<=b; k++) x = mulsi(k,x);
    3346         287 :     return gerepileuptoint(av, x);
    3347             :   }
    3348           0 :   lx = 1; x = cgetg(2 + n/2, t_VEC);
    3349           0 :   N = b + a;
    3350           0 :   for (k = a;; k++)
    3351             :   {
    3352           0 :     l = N - k; if (l <= k) break;
    3353           0 :     gel(x,lx++) = mulss(k,l);
    3354             :   }
    3355           0 :   if (l == k) gel(x,lx++) = stoi(k);
    3356           0 :   setlg(x, lx);
    3357           0 :   return gerepileuptoint(av, ZV_prod(x));
    3358             : }
    3359             : 
    3360             : GEN
    3361         105 : mpprimorial(long n)
    3362             : {
    3363         105 :   pari_sp av = avma;
    3364         105 :   if (n <= 12) switch(n)
    3365             :   {
    3366          14 :     case 0: case 1: return gen_1;
    3367           7 :     case 2: return gen_2;
    3368          14 :     case 3: case 4: return utoipos(6);
    3369          14 :     case 5: case 6: return utoipos(30);
    3370          28 :     case 7: case 8: case 9: case 10: return utoipos(210);
    3371          14 :     case 11: case 12: return utoipos(2310);
    3372           7 :     default: pari_err_DOMAIN("primorial", "argument","<",gen_0,stoi(n));
    3373             :   }
    3374           7 :   return gerepileuptoint(av, zv_prod_Z(primes_upto_zv(n)));
    3375             : }
    3376             : 
    3377             : GEN
    3378      496582 : mpfact(long n)
    3379             : {
    3380      496582 :   pari_sp av = avma;
    3381             :   GEN a, v;
    3382             :   long k;
    3383      496582 :   if (n <= 12) switch(n)
    3384             :   {
    3385      428654 :     case 0: case 1: return gen_1;
    3386       24338 :     case 2: return gen_2;
    3387        3388 :     case 3: return utoipos(6);
    3388        4145 :     case 4: return utoipos(24);
    3389        2887 :     case 5: return utoipos(120);
    3390        2556 :     case 6: return utoipos(720);
    3391        2448 :     case 7: return utoipos(5040);
    3392        2437 :     case 8: return utoipos(40320);
    3393        2458 :     case 9: return utoipos(362880);
    3394        2694 :     case 10:return utoipos(3628800);
    3395        1409 :     case 11:return utoipos(39916800);
    3396         577 :     case 12:return utoipos(479001600);
    3397           0 :     default: pari_err_DOMAIN("factorial", "argument","<",gen_0,stoi(n));
    3398             :   }
    3399       18591 :   v = cgetg(expu(n) + 2, t_VEC);
    3400       18578 :   for (k = 1;; k++)
    3401       86842 :   {
    3402      105420 :     long m = n >> (k-1), l;
    3403      105420 :     if (m <= 2) break;
    3404       86833 :     l = (1 + (n >> k)) | 1;
    3405             :     /* product of odd numbers in ]n / 2^k, n / 2^(k-1)] */
    3406       86833 :     a = mulu_interval_step(l, m, 2);
    3407       86793 :     gel(v,k) = k == 1? a: powiu(a, k);
    3408             :   }
    3409       86874 :   a = gel(v,--k); while (--k) a = mulii(a, gel(v,k));
    3410       18588 :   a = shifti(a, factorial_lval(n, 2));
    3411       18585 :   return gerepileuptoint(av, a);
    3412             : }
    3413             : 
    3414             : ulong
    3415       56831 : factorial_Fl(long n, ulong p)
    3416             : {
    3417             :   long k;
    3418             :   ulong v;
    3419       56831 :   if (p <= (ulong)n) return 0;
    3420       56831 :   v = Fl_powu(2, factorial_lval(n, 2), p);
    3421       56891 :   for (k = 1;; k++)
    3422      142604 :   {
    3423      199495 :     long m = n >> (k-1), l, i;
    3424      199495 :     ulong a = 1;
    3425      199495 :     if (m <= 2) break;
    3426      142610 :     l = (1 + (n >> k)) | 1;
    3427             :     /* product of odd numbers in ]n / 2^k, 2 / 2^(k-1)] */
    3428      779310 :     for (i=l; i<=m; i+=2)
    3429      636700 :       a = Fl_mul(a, i, p);
    3430      142610 :     v = Fl_mul(v, k == 1? a: Fl_powu(a, k, p), p);
    3431             :   }
    3432       56885 :   return v;
    3433             : }
    3434             : 
    3435             : GEN
    3436         382 : factorial_Fp(long n, GEN p)
    3437             : {
    3438         382 :   pari_sp av = avma;
    3439             :   long k;
    3440         382 :   GEN v = Fp_powu(gen_2, factorial_lval(n, 2), p);
    3441         382 :   for (k = 1;; k++)
    3442        1240 :   {
    3443        1622 :     long m = n >> (k-1), l, i;
    3444        1622 :     GEN a = gen_1;
    3445        1622 :     if (m <= 2) break;
    3446        1240 :     l = (1 + (n >> k)) | 1;
    3447             :     /* product of odd numbers in ]n / 2^k, 2 / 2^(k-1)] */
    3448        7570 :     for (i=l; i<=m; i+=2)
    3449        6330 :       a = Fp_mulu(a, i, p);
    3450        1240 :     v = Fp_mul(v, k == 1? a: Fp_powu(a, k, p), p);
    3451        1240 :     v = gerepileuptoint(av, v);
    3452             :   }
    3453         382 :   return v;
    3454             : }
    3455             : 
    3456             : /*******************************************************************/
    3457             : /**                      LUCAS & FIBONACCI                        **/
    3458             : /*******************************************************************/
    3459             : static void
    3460          56 : lucas(ulong n, GEN *a, GEN *b)
    3461             : {
    3462             :   GEN z, t, zt;
    3463          56 :   if (!n) { *a = gen_2; *b = gen_1; return; }
    3464          49 :   lucas(n >> 1, &z, &t); zt = mulii(z, t);
    3465          49 :   switch(n & 3) {
    3466          14 :     case  0: *a = subiu(sqri(z),2); *b = subiu(zt,1); break;
    3467          14 :     case  1: *a = subiu(zt,1);      *b = addiu(sqri(t),2); break;
    3468           7 :     case  2: *a = addiu(sqri(z),2); *b = addiu(zt,1); break;
    3469          14 :     case  3: *a = addiu(zt,1);      *b = subiu(sqri(t),2);
    3470             :   }
    3471             : }
    3472             : 
    3473             : GEN
    3474           7 : fibo(long n)
    3475             : {
    3476           7 :   pari_sp av = avma;
    3477             :   GEN a, b;
    3478           7 :   if (!n) return gen_0;
    3479           7 :   lucas((ulong)(labs(n)-1), &a, &b);
    3480           7 :   a = diviuexact(addii(shifti(a,1),b), 5);
    3481           7 :   if (n < 0 && !odd(n)) setsigne(a, -1);
    3482           7 :   return gerepileuptoint(av, a);
    3483             : }
    3484             : 
    3485             : /*******************************************************************/
    3486             : /*                      CONTINUED FRACTIONS                        */
    3487             : /*******************************************************************/
    3488             : static GEN
    3489     3136994 : icopy_lg(GEN x, long l)
    3490             : {
    3491     3136994 :   long lx = lgefint(x);
    3492             :   GEN y;
    3493             : 
    3494     3136994 :   if (lx >= l) return icopy(x);
    3495          49 :   y = cgeti(l); affii(x, y); return y;
    3496             : }
    3497             : 
    3498             : /* continued fraction of a/b. If y != NULL, stop when partial quotients
    3499             :  * differ from y */
    3500             : static GEN
    3501     3137344 : Qsfcont(GEN a, GEN b, GEN y, ulong k)
    3502             : {
    3503             :   GEN  z, c;
    3504     3137344 :   ulong i, l, ly = lgefint(b);
    3505             : 
    3506             :   /* times 1 / log2( (1+sqrt(5)) / 2 )  */
    3507     3137344 :   l = (ulong)(3 + bit_accuracy_mul(ly, 1.44042009041256));
    3508     3137344 :   if (k > 0 && k+1 > 0 && l > k+1) l = k+1; /* beware overflow */
    3509     3137344 :   if (l > LGBITS) l = LGBITS;
    3510             : 
    3511     3137344 :   z = cgetg(l,t_VEC);
    3512     3137344 :   l--;
    3513     3137344 :   if (y) {
    3514         350 :     pari_sp av = avma;
    3515         350 :     if (l >= (ulong)lg(y)) l = lg(y)-1;
    3516       25209 :     for (i = 1; i <= l; i++)
    3517             :     {
    3518       24985 :       GEN q = gel(y,i);
    3519       24985 :       gel(z,i) = q;
    3520       24985 :       c = b; if (!gequal1(q)) c = mulii(q, b);
    3521       24985 :       c = subii(a, c);
    3522       24985 :       if (signe(c) < 0)
    3523             :       { /* partial quotient too large */
    3524          96 :         c = addii(c, b);
    3525          96 :         if (signe(c) >= 0) i++; /* by 1 */
    3526          96 :         break;
    3527             :       }
    3528       24889 :       if (cmpii(c, b) >= 0)
    3529             :       { /* partial quotient too small */
    3530          30 :         c = subii(c, b);
    3531          30 :         if (cmpii(c, b) < 0) {
    3532             :           /* by 1. If next quotient is 1 in y, add 1 */
    3533          12 :           if (i < l && equali1(gel(y,i+1))) gel(z,i) = addiu(q,1);
    3534          12 :           i++;
    3535             :         }
    3536          30 :         break;
    3537             :       }
    3538       24859 :       if ((i & 0xff) == 0) gerepileall(av, 2, &b, &c);
    3539       24859 :       a = b; b = c;
    3540             :     }
    3541             :   } else {
    3542     3136994 :     a = icopy_lg(a, ly);
    3543     3136994 :     b = icopy(b);
    3544    24524282 :     for (i = 1; i <= l; i++)
    3545             :     {
    3546    24523964 :       gel(z,i) = truedvmdii(a,b,&c);
    3547    24523964 :       if (c == gen_0) { i++; break; }
    3548    21387288 :       affii(c, a); cgiv(c); c = a;
    3549    21387288 :       a = b; b = c;
    3550             :     }
    3551             :   }
    3552     3137344 :   i--;
    3553     3137344 :   if (i > 1 && gequal1(gel(z,i)))
    3554             :   {
    3555         101 :     cgiv(gel(z,i)); --i;
    3556         101 :     gel(z,i) = addui(1, gel(z,i)); /* unclean: leave old z[i] on stack */
    3557             :   }
    3558     3137344 :   setlg(z,i+1); return z;
    3559             : }
    3560             : 
    3561             : static GEN
    3562           0 : sersfcont(GEN a, GEN b, long k)
    3563             : {
    3564           0 :   long i, l = typ(a) == t_POL? lg(a): 3;
    3565             :   GEN y, c;
    3566           0 :   if (lg(b) > l) l = lg(b);
    3567           0 :   if (k > 0 && l > k+1) l = k+1;
    3568           0 :   y = cgetg(l,t_VEC);
    3569           0 :   for (i=1; i<l; i++)
    3570             :   {
    3571           0 :     gel(y,i) = poldivrem(a,b,&c);
    3572           0 :     if (gequal0(c)) { i++; break; }
    3573           0 :     a = b; b = c;
    3574             :   }
    3575           0 :   setlg(y, i); return y;
    3576             : }
    3577             : 
    3578             : GEN
    3579     3142307 : gboundcf(GEN x, long k)
    3580             : {
    3581             :   pari_sp av;
    3582     3142307 :   long tx = typ(x), e;
    3583             :   GEN y, a, b, c;
    3584             : 
    3585     3142307 :   if (k < 0) pari_err_DOMAIN("gboundcf","nmax","<",gen_0,stoi(k));
    3586     3142300 :   if (is_scalar_t(tx))
    3587             :   {
    3588     3142300 :     if (gequal0(x)) return mkvec(gen_0);
    3589     3142181 :     switch(tx)
    3590             :     {
    3591        5180 :       case t_INT: return mkveccopy(x);
    3592         357 :       case t_REAL:
    3593         357 :         av = avma;
    3594         357 :         c = mantissa_real(x,&e);
    3595         357 :         if (e < 0) pari_err_PREC("gboundcf");
    3596         350 :         y = int2n(e);
    3597         350 :         a = Qsfcont(c,y, NULL, k);
    3598         350 :         b = addsi(signe(x), c);
    3599         350 :         return gerepilecopy(av, Qsfcont(b,y, a, k));
    3600             : 
    3601     3136644 :       case t_FRAC:
    3602     3136644 :         av = avma;
    3603     3136644 :         return gerepileupto(av, Qsfcont(gel(x,1),gel(x,2), NULL, k));
    3604             :     }
    3605           0 :     pari_err_TYPE("gboundcf",x);
    3606             :   }
    3607             : 
    3608           0 :   switch(tx)
    3609             :   {
    3610           0 :     case t_POL: return mkveccopy(x);
    3611           0 :     case t_SER:
    3612           0 :       av = avma;
    3613           0 :       return gerepileupto(av, gboundcf(ser2rfrac_i(x), k));
    3614           0 :     case t_RFRAC:
    3615           0 :       av = avma;
    3616           0 :       return gerepilecopy(av, sersfcont(gel(x,1), gel(x,2), k));
    3617             :   }
    3618           0 :   pari_err_TYPE("gboundcf",x);
    3619             :   return NULL; /* LCOV_EXCL_LINE */
    3620             : }
    3621             : 
    3622             : static GEN
    3623          14 : sfcont2(GEN b, GEN x, long k)
    3624             : {
    3625          14 :   pari_sp av = avma;
    3626          14 :   long lb = lg(b), tx = typ(x), i;
    3627             :   GEN y,p1;
    3628             : 
    3629          14 :   if (k)
    3630             :   {
    3631           7 :     if (k >= lb) pari_err_DIM("contfrac [too few denominators]");
    3632           0 :     lb = k+1;
    3633             :   }
    3634           7 :   y = cgetg(lb,t_VEC);
    3635           7 :   if (lb==1) return y;
    3636           7 :   if (is_scalar_t(tx))
    3637             :   {
    3638           7 :     if (!is_intreal_t(tx) && tx != t_FRAC) pari_err_TYPE("sfcont2",x);
    3639             :   }
    3640           0 :   else if (tx == t_SER) x = ser2rfrac_i(x);
    3641             : 
    3642           7 :   if (!gequal1(gel(b,1))) x = gmul(gel(b,1),x);
    3643           7 :   for (i = 1;;)
    3644             :   {
    3645          35 :     if (tx == t_REAL)
    3646             :     {
    3647          35 :       long e = expo(x);
    3648          35 :       if (e > 0 && nbits2prec(e+1) > realprec(x)) break;
    3649          35 :       gel(y,i) = floorr(x);
    3650          35 :       p1 = subri(x, gel(y,i));
    3651             :     }
    3652             :     else
    3653             :     {
    3654           0 :       gel(y,i) = gfloor(x);
    3655           0 :       p1 = gsub(x, gel(y,i));
    3656             :     }
    3657          35 :     if (++i >= lb) break;
    3658          28 :     if (gequal0(p1)) break;
    3659          28 :     x = gdiv(gel(b,i),p1);
    3660             :   }
    3661           7 :   setlg(y,i);
    3662           7 :   return gerepilecopy(av,y);
    3663             : }
    3664             : 
    3665             : GEN
    3666         126 : gcf(GEN x) { return gboundcf(x,0); }
    3667             : GEN
    3668           0 : gcf2(GEN b, GEN x) { return contfrac0(x,b,0); }
    3669             : GEN
    3670          49 : contfrac0(GEN x, GEN b, long nmax)
    3671             : {
    3672             :   long tb;
    3673             : 
    3674          49 :   if (!b) return gboundcf(x,nmax);
    3675          28 :   tb = typ(b);
    3676          28 :   if (tb == t_INT) return gboundcf(x,itos(b));
    3677          21 :   if (! is_vec_t(tb)) pari_err_TYPE("contfrac0",b);
    3678          21 :   if (nmax < 0) pari_err_DOMAIN("contfrac","nmax","<",gen_0,stoi(nmax));
    3679          14 :   return sfcont2(b,x,nmax);
    3680             : }
    3681             : 
    3682             : GEN
    3683         266 : contfracpnqn(GEN x, long n)
    3684             : {
    3685         266 :   pari_sp av = avma;
    3686         266 :   long i, lx = lg(x);
    3687             :   GEN M,A,B, p0,p1, q0,q1;
    3688             : 
    3689         266 :   if (lx == 1)
    3690             :   {
    3691          28 :     if (! is_matvec_t(typ(x))) pari_err_TYPE("pnqn",x);
    3692          21 :     if (n >= 0) return cgetg(1,t_MAT);
    3693           7 :     return matid(2);
    3694             :   }
    3695         238 :   switch(typ(x))
    3696             :   {
    3697         196 :     case t_VEC: case t_COL: A = x; B = NULL; break;
    3698          42 :     case t_MAT:
    3699          42 :       switch(lgcols(x))
    3700             :       {
    3701           0 :         case 2: A = row(x,1); B = NULL; break;
    3702          35 :         case 3: A = row(x,2); B = row(x,1); break;
    3703           7 :         default: pari_err_DIM("pnqn [ nbrows != 1,2 ]");
    3704             :                  return NULL; /*LCOV_EXCL_LINE*/
    3705             :       }
    3706          35 :       break;
    3707           0 :     default: pari_err_TYPE("pnqn",x);
    3708             :       return NULL; /*LCOV_EXCL_LINE*/
    3709             :   }
    3710         231 :   p1 = gel(A,1);
    3711         231 :   q1 = B? gel(B,1): gen_1; /* p[0], q[0] */
    3712         231 :   if (n >= 0)
    3713             :   {
    3714         196 :     lx = minss(lx, n+2);
    3715         196 :     if (lx == 2) return gerepilecopy(av, mkmat(mkcol2(p1,q1)));
    3716             :   }
    3717          35 :   else if (lx == 2)
    3718           7 :     return gerepilecopy(av, mkmat2(mkcol2(p1,q1), mkcol2(gen_1,gen_0)));
    3719             :   /* lx >= 3 */
    3720         119 :   p0 = gen_1;
    3721         119 :   q0 = gen_0; /* p[-1], q[-1] */
    3722         119 :   M = cgetg(lx, t_MAT);
    3723         119 :   gel(M,1) = mkcol2(p1,q1);
    3724         399 :   for (i=2; i<lx; i++)
    3725             :   {
    3726         280 :     GEN a = gel(A,i), p2,q2;
    3727         280 :     if (B) {
    3728          84 :       GEN b = gel(B,i);
    3729          84 :       p0 = gmul(b,p0);
    3730          84 :       q0 = gmul(b,q0);
    3731             :     }
    3732         280 :     p2 = gadd(gmul(a,p1),p0); p0=p1; p1=p2;
    3733         280 :     q2 = gadd(gmul(a,q1),q0); q0=q1; q1=q2;
    3734         280 :     gel(M,i) = mkcol2(p1,q1);
    3735             :   }
    3736         119 :   if (n < 0) M = mkmat2(gel(M,lx-1), gel(M,lx-2));
    3737         119 :   return gerepilecopy(av, M);
    3738             : }
    3739             : GEN
    3740           0 : pnqn(GEN x) { return contfracpnqn(x,-1); }
    3741             : /* x = [a0, ..., an] from gboundcf, n >= 0;
    3742             :  * return [[p0, ..., pn], [q0,...,qn]] */
    3743             : GEN
    3744      894782 : ZV_allpnqn(GEN x)
    3745             : {
    3746      894782 :   long i, lx = lg(x);
    3747      894782 :   GEN p0, p1, q0, q1, p2, q2, P,Q, v = cgetg(3,t_VEC);
    3748             : 
    3749      894782 :   gel(v,1) = P = cgetg(lx, t_VEC);
    3750      894782 :   gel(v,2) = Q = cgetg(lx, t_VEC);
    3751      894782 :   p0 = gen_1; q0 = gen_0;
    3752      894782 :   gel(P, 1) = p1 = gel(x,1); gel(Q, 1) = q1 = gen_1;
    3753     3106138 :   for (i=2; i<lx; i++)
    3754             :   {
    3755     2211356 :     GEN a = gel(x,i);
    3756     2211356 :     gel(P, i) = p2 = addmulii(p0, a, p1); p0 = p1; p1 = p2;
    3757     2211356 :     gel(Q, i) = q2 = addmulii(q0, a, q1); q0 = q1; q1 = q2;
    3758             :   }
    3759      894782 :   return v;
    3760             : }
    3761             : 
    3762             : /* write Mod(x,N) as a/b, gcd(a,b) = 1, b <= B (no condition if B = NULL) */
    3763             : static GEN
    3764          42 : mod_to_frac(GEN x, GEN N, GEN B)
    3765             : {
    3766             :   GEN a, b, A;
    3767          42 :   if (B) A = divii(shifti(N, -1), B);
    3768             :   else
    3769             :   {
    3770          14 :     A = sqrti(shifti(N, -1));
    3771          14 :     B = A;
    3772             :   }
    3773          42 :   if (!Fp_ratlift(x, N, A,B,&a,&b) || !equali1( gcdii(a,b) )) return NULL;
    3774          28 :   return equali1(b)? a: mkfrac(a,b);
    3775             : }
    3776             : 
    3777             : static GEN
    3778         112 : mod_to_rfrac(GEN x, GEN N, long B)
    3779             : {
    3780             :   GEN a, b;
    3781         112 :   long A, d = degpol(N);
    3782         112 :   if (B >= 0) A = d-1 - B;
    3783             :   else
    3784             :   {
    3785          42 :     B = d >> 1;
    3786          42 :     A = odd(d)? B : B-1;
    3787             :   }
    3788         112 :   if (varn(N) != varn(x)) x = scalarpol(x, varn(N));
    3789         112 :   if (!RgXQ_ratlift(x, N, A, B, &a,&b) || degpol(RgX_gcd(a,b)) > 0) return NULL;
    3790          91 :   return gdiv(a,b);
    3791             : }
    3792             : 
    3793             : /* k > 0 t_INT, x a t_FRAC, returns the convergent a/b
    3794             :  * of the continued fraction of x with b <= k maximal */
    3795             : static GEN
    3796           7 : bestappr_frac(GEN x, GEN k)
    3797             : {
    3798             :   pari_sp av;
    3799             :   GEN p0, p1, p, q0, q1, q, a, y;
    3800             : 
    3801           7 :   if (cmpii(gel(x,2),k) <= 0) return gcopy(x);
    3802           0 :   av = avma; y = x;
    3803           0 :   p1 = gen_1; p0 = truedvmdii(gel(x,1), gel(x,2), &a); /* = floor(x) */
    3804           0 :   q1 = gen_0; q0 = gen_1;
    3805           0 :   x = mkfrac(a, gel(x,2)); /* = frac(x); now 0<= x < 1 */
    3806             :   for(;;)
    3807             :   {
    3808           0 :     x = ginv(x); /* > 1 */
    3809           0 :     a = typ(x)==t_INT? x: divii(gel(x,1), gel(x,2));
    3810           0 :     if (cmpii(a,k) > 0)
    3811             :     { /* next partial quotient will overflow limits */
    3812             :       GEN n, d;
    3813           0 :       a = divii(subii(k, q1), q0);
    3814           0 :       p = addii(mulii(a,p0), p1); p1=p0; p0=p;
    3815           0 :       q = addii(mulii(a,q0), q1); q1=q0; q0=q;
    3816             :       /* compare |y-p0/q0|, |y-p1/q1| */
    3817           0 :       n = gel(y,1);
    3818           0 :       d = gel(y,2);
    3819           0 :       if (abscmpii(mulii(q1, subii(mulii(q0,n), mulii(d,p0))),
    3820             :                    mulii(q0, subii(mulii(q1,n), mulii(d,p1)))) < 0)
    3821           0 :                    { p1 = p0; q1 = q0; }
    3822           0 :       break;
    3823             :     }
    3824           0 :     p = addii(mulii(a,p0), p1); p1=p0; p0=p;
    3825           0 :     q = addii(mulii(a,q0), q1); q1=q0; q0=q;
    3826             : 
    3827           0 :     if (cmpii(q0,k) > 0) break;
    3828           0 :     x = gsub(x,a); /* 0 <= x < 1 */
    3829           0 :     if (typ(x) == t_INT) { p1 = p0; q1 = q0; break; } /* x = 0 */
    3830             : 
    3831             :   }
    3832           0 :   return gerepileupto(av, gdiv(p1,q1));
    3833             : }
    3834             : /* k > 0 t_INT, x != 0 a t_REAL, returns the convergent a/b
    3835             :  * of the continued fraction of x with b <= k maximal */
    3836             : static GEN
    3837     1267968 : bestappr_real(GEN x, GEN k)
    3838             : {
    3839     1267968 :   pari_sp av = avma;
    3840     1267968 :   GEN kr, p0, p1, p, q0, q1, q, a, y = x;
    3841             : 
    3842     1267968 :   p1 = gen_1; a = p0 = floorr(x);
    3843     1267880 :   q1 = gen_0; q0 = gen_1;
    3844     1267880 :   x = subri(x,a); /* 0 <= x < 1 */
    3845     1267901 :   if (!signe(x)) { cgiv(x); return a; }
    3846     1151639 :   kr = itor(k, realprec(x));
    3847             :   for(;;)
    3848     1277387 :   {
    3849             :     long d;
    3850     2429054 :     x = invr(x); /* > 1 */
    3851     2428813 :     if (cmprr(x,kr) > 0)
    3852             :     { /* next partial quotient will overflow limits */
    3853     1129639 :       a = divii(subii(k, q1), q0);
    3854     1129637 :       p = addii(mulii(a,p0), p1); p1=p0; p0=p;
    3855     1129670 :       q = addii(mulii(a,q0), q1); q1=q0; q0=q;
    3856             :       /* compare |y-p0/q0|, |y-p1/q1| */
    3857     1129644 :       if (abscmprr(mulir(q1, subri(mulir(q0,y), p0)),
    3858             :                    mulir(q0, subri(mulir(q1,y), p1))) < 0)
    3859      125959 :                    { p1 = p0; q1 = q0; }
    3860     1129704 :       break;
    3861             :     }
    3862     1299278 :     d = nbits2prec(expo(x) + 1);
    3863     1299278 :     if (d > realprec(x)) { p1 = p0; q1 = q0; break; } /* original x was ~ 0 */
    3864             : 
    3865     1299088 :     a = truncr(x); /* truncr(x) will NOT raise e_PREC */
    3866     1299049 :     p = addii(mulii(a,p0), p1); p1=p0; p0=p;
    3867     1299071 :     q = addii(mulii(a,q0), q1); q1=q0; q0=q;
    3868             : 
    3869     1299074 :     if (cmpii(q0,k) > 0) break;
    3870     1292819 :     x = subri(x,a); /* 0 <= x < 1 */
    3871     1292835 :     if (!signe(x)) { p1 = p0; q1 = q0; break; }
    3872             :   }
    3873     1151597 :   if (signe(q1) < 0) { togglesign_safe(&p1); togglesign_safe(&q1); }
    3874     1151597 :   return gerepilecopy(av, equali1(q1)? p1: mkfrac(p1,q1));
    3875             : }
    3876             : 
    3877             : /* k t_INT or NULL */
    3878             : static GEN
    3879     2271587 : bestappr_Q(GEN x, GEN k)
    3880             : {
    3881     2271587 :   long lx, tx = typ(x), i;
    3882             :   GEN a, y;
    3883             : 
    3884     2271587 :   switch(tx)
    3885             :   {
    3886         154 :     case t_INT: return icopy(x);
    3887           7 :     case t_FRAC: return k? bestappr_frac(x, k): gcopy(x);
    3888     1521956 :     case t_REAL:
    3889     1521956 :       if (!signe(x)) return gen_0;
    3890             :       /* i <= e iff nbits2lg(e+1) > lg(x) iff floorr(x) fails */
    3891     1267955 :       i = bit_prec(x); if (i <= expo(x)) return NULL;
    3892     1267970 :       return bestappr_real(x, k? k: int2n(i));
    3893             : 
    3894          28 :     case t_INTMOD: {
    3895          28 :       pari_sp av = avma;
    3896          28 :       a = mod_to_frac(gel(x,2), gel(x,1), k); if (!a) return NULL;
    3897          21 :       return gerepilecopy(av, a);
    3898             :     }
    3899          14 :     case t_PADIC: {
    3900          14 :       pari_sp av = avma;
    3901          14 :       long v = valp(x);
    3902          14 :       a = mod_to_frac(padic_u(x), padic_pd(x), k); if (!a) return NULL;
    3903           7 :       if (v) a = gmul(a, powis(padic_p(x), v));
    3904           7 :       return gerepilecopy(av, a);
    3905             :     }
    3906             : 
    3907        5453 :     case t_COMPLEX: {
    3908        5453 :       pari_sp av = avma;
    3909        5453 :       y = cgetg(3, t_COMPLEX);
    3910        5453 :       gel(y,2) = bestappr(gel(x,2), k);
    3911        5453 :       gel(y,1) = bestappr(gel(x,1), k);
    3912        5453 :       if (gequal0(gel(y,2))) return gerepileupto(av, gel(y,1));
    3913          91 :       return y;
    3914             :     }
    3915           0 :     case t_SER:
    3916           0 :       if (ser_isexactzero(x)) return gcopy(x);
    3917             :       /* fall through */
    3918             :     case t_POLMOD: case t_POL: case t_RFRAC:
    3919             :     case t_VEC: case t_COL: case t_MAT:
    3920      743975 :       y = cgetg_copy(x, &lx);
    3921      744023 :       for(i = 1; i < lontyp[tx]; i++) y[i] = x[i];
    3922     2910483 :       for (; i < lx; i++)
    3923             :       {
    3924     2166503 :         a = bestappr_Q(gel(x,i),k); if (!a) return NULL;
    3925     2166474 :         gel(y,i) = a;
    3926             :       }
    3927      743980 :       if (tx == t_POL) return normalizepol(y);
    3928      743966 :       if (tx == t_SER) return normalizeser(y);
    3929      743966 :       return y;
    3930             :   }
    3931           0 :   pari_err_TYPE("bestappr_Q",x);
    3932             :   return NULL; /* LCOV_EXCL_LINE */
    3933             : }
    3934             : 
    3935             : static GEN
    3936          98 : bestappr_ser(GEN x, long B)
    3937             : {
    3938          98 :   long dN, v = valser(x), lx = lg(x);
    3939             :   GEN t;
    3940          98 :   x = normalizepol(ser2pol_i(x, lx));
    3941          98 :   dN = lx-2;
    3942          98 :   if (v > 0)
    3943             :   {
    3944          21 :     x = RgX_shift_shallow(x, v);
    3945          21 :     dN += v;
    3946             :   }
    3947          77 :   else if (v < 0)
    3948             :   {
    3949          14 :     if (B >= 0) B = maxss(B+v, 0);
    3950             :   }
    3951          98 :   t = mod_to_rfrac(x, pol_xn(dN, varn(x)), B);
    3952          98 :   if (!t) return NULL;
    3953          77 :   if (v < 0)
    3954             :   {
    3955             :     GEN a, b;
    3956             :     long vx;
    3957          14 :     if (typ(t) == t_POL) return RgX_mulXn(t, v);
    3958             :     /* t_RFRAC */
    3959          14 :     vx = varn(x);
    3960          14 :     a = gel(t,1);
    3961          14 :     b = gel(t,2);
    3962          14 :     v -= RgX_valrem(b, &b);
    3963          14 :     if (typ(a) == t_POL && varn(a) == vx) v += RgX_valrem(a, &a);
    3964          14 :     if (v < 0) b = RgX_shift_shallow(b, -v);
    3965           0 :     else if (v > 0) {
    3966           0 :       if (typ(a) != t_POL || varn(a) != vx) a = scalarpol_shallow(a, vx);
    3967           0 :       a = RgX_shift_shallow(a, v);
    3968             :     }
    3969          14 :     t = mkrfraccopy(a, b);
    3970             :   }
    3971          77 :   return t;
    3972             : }
    3973             : static GEN
    3974          42 : gc_empty(pari_sp av) { set_avma(av); return cgetg(1, t_VEC); }
    3975             : static GEN
    3976         112 : _gc_upto(pari_sp av, GEN x) { return x? gerepileupto(av, x): NULL; }
    3977             : 
    3978             : static GEN bestappr_RgX(GEN x, long B);
    3979             : /* B >= 0 or < 0 [omit condition on B].
    3980             :  * Look for coprime t_POL a,b, deg(b)<=B, such that a/b ~ x */
    3981             : static GEN
    3982         119 : bestappr_RgX(GEN x, long B)
    3983             : {
    3984             :   pari_sp av;
    3985         119 :   switch(typ(x))
    3986             :   {
    3987           0 :     case t_INT: case t_REAL: case t_INTMOD: case t_FRAC: case t_FFELT:
    3988             :     case t_COMPLEX: case t_PADIC: case t_QUAD: case t_POL:
    3989           0 :       return gcopy(x);
    3990          14 :     case t_RFRAC:
    3991          14 :       if (B < 0 || degpol(gel(x,2)) <= B) return gcopy(x);
    3992           7 :       av = avma; return _gc_upto(av, bestappr_ser(rfrac_to_ser_i(x, 2*B+1), B));
    3993          14 :     case t_POLMOD:
    3994          14 :       av = avma; return _gc_upto(av, mod_to_rfrac(gel(x,2), gel(x,1), B));
    3995          91 :     case t_SER:
    3996          91 :       av = avma; return _gc_upto(av, bestappr_ser(x, B));
    3997           0 :     case t_VEC: case t_COL: case t_MAT: {
    3998             :       long i, lx;
    3999           0 :       GEN y = cgetg_copy(x, &lx);
    4000           0 :       for (i = 1; i < lx; i++)
    4001             :       {
    4002           0 :         GEN t = bestappr_RgX(gel(x,i),B); if (!t) return NULL;
    4003           0 :         gel(y,i) = t;
    4004             :       }
    4005           0 :       return y;
    4006             :     }
    4007             :   }
    4008           0 :   pari_err_TYPE("bestappr_RgX",x);
    4009             :   return NULL; /* LCOV_EXCL_LINE */
    4010             : }
    4011             : 
    4012             : /* allow k = NULL: maximal accuracy */
    4013             : GEN
    4014      105075 : bestappr(GEN x, GEN k)
    4015             : {
    4016      105075 :   pari_sp av = avma;
    4017      105075 :   if (k) { /* replace by floor(k) */
    4018      104753 :     switch(typ(k))
    4019             :     {
    4020       33026 :       case t_INT:
    4021       33026 :         break;
    4022       71727 :       case t_REAL: case t_FRAC:
    4023       71727 :         k = floor_safe(k); /* left on stack for efficiency */
    4024       71726 :         if (!signe(k)) k = gen_1;
    4025       71726 :         break;
    4026           0 :       default:
    4027           0 :         pari_err_TYPE("bestappr [bound type]", k);
    4028           0 :         break;
    4029             :     }
    4030             :   }
    4031      105074 :   x = bestappr_Q(x, k);
    4032      105074 :   return x? x: gc_empty(av);
    4033             : }
    4034             : GEN
    4035         119 : bestapprPade(GEN x, long B)
    4036             : {
    4037         119 :   pari_sp av = avma;
    4038         119 :   GEN t = bestappr_RgX(x, B);
    4039         119 :   return t? t: gc_empty(av);
    4040             : }
    4041             : 
    4042             : static GEN
    4043          49 : serPade(GEN S, long p, long q)
    4044             : {
    4045          49 :   pari_sp av = avma;
    4046          49 :   long va, v, t = typ(S);
    4047          49 :   if (t!=t_SER && t!=t_POL && t!=t_RFRAC) pari_err_TYPE("bestapprPade", S);
    4048          49 :   va = gvar(S); v = gvaluation(S, pol_x(va));
    4049          49 :   if (p < 0) pari_err_DOMAIN("bestapprPade", "p", "<", gen_0, stoi(p));
    4050          49 :   if (q < 0) pari_err_DOMAIN("bestapprPade", "q", "<", gen_0, stoi(q));
    4051          49 :   if (v == LONG_MAX) return gc_empty(av);
    4052          42 :   S = gadd(S, zeroser(va, p + q + 1 + v));
    4053          42 :   return gerepileupto(av, bestapprPade(S, q));
    4054             : }
    4055             : 
    4056             : GEN
    4057         126 : bestapprPade0(GEN x, long p, long q)
    4058             : {
    4059          77 :   return (p >= 0 && q >= 0)? serPade(x, p, q)
    4060         203 :                            : bestapprPade(x, p >= 0? p: q);
    4061             : }

Generated by: LCOV version 1.16