Code coverage tests

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

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

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

LCOV - code coverage report
Current view: top level - basemath - arith1.c (source / functions) Hit Total Coverage
Test: PARI/GP v2.16.1 lcov report (development 28574-61f195acfe) Lines: 2016 2193 91.9 %
Date: 2023-06-08 07:47:28 Functions: 209 220 95.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      464816 : 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      464816 : 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        4895 : is_gener_expo(GEN p, GEN L0)
      39             : {
      40        4895 :   GEN L, q = shifti(p,-1);
      41             :   long i, l;
      42        4895 :   if (L0) {
      43        3120 :     l = lg(L0);
      44        3120 :     L = cgetg(l, t_VEC);
      45             :   } else {
      46        1775 :     L0 = L = odd_prime_divisors(q);
      47        1775 :     l = lg(L);
      48             :   }
      49       14171 :   for (i=1; i<l; i++) gel(L,i) = diviiexact(q, gel(L0,i));
      50        4895 :   return L;
      51             : }
      52             : static GEN
      53      532669 : u_is_gener_expo(ulong p, GEN L0)
      54             : {
      55      532669 :   const ulong q = p >> 1;
      56             :   long i;
      57             :   GEN L;
      58      532669 :   if (!L0) L0 = u_odd_prime_divisors(q);
      59      532671 :   L = cgetg_copy(L0,&i);
      60     1154949 :   while (--i) L[i] = q / uel(L0,i);
      61      532670 :   return L;
      62             : }
      63             : 
      64             : int
      65     1629492 : is_gener_Fl(ulong x, ulong p, ulong p_1, GEN L)
      66             : {
      67             :   long i;
      68     1629492 :   if (krouu(x, p) >= 0) return 0;
      69     1379928 :   for (i=lg(L)-1; i; i--)
      70             :   {
      71      842462 :     ulong t = Fl_powu(x, uel(L,i), p);
      72      842460 :     if (t == p_1 || t == 1) return 0;
      73             :   }
      74      537466 :   return 1;
      75             : }
      76             : /* assume p prime */
      77             : ulong
      78     1054026 : pgener_Fl_local(ulong p, GEN L0)
      79             : {
      80     1054026 :   const pari_sp av = avma;
      81     1054026 :   const ulong p_1 = p-1;
      82             :   long x;
      83             :   GEN L;
      84     1054026 :   if (p <= 19) switch(p)
      85             :   { /* quick trivial cases */
      86          28 :     case 2:  return 1;
      87      116176 :     case 7:
      88      116176 :     case 17: return 3;
      89      405199 :     default: return 2;
      90             :   }
      91      532623 :   L = u_is_gener_expo(p,L0);
      92     1621675 :   for (x = 2;; x++)
      93     1621675 :     if (is_gener_Fl(x,p,p_1,L)) return gc_ulong(av, x);
      94             : }
      95             : ulong
      96      579292 : 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       13755 : is_gener_Fp(GEN x, GEN p, GEN p_1, GEN L)
     102             : {
     103       13755 :   long i, t = lgefint(x)==3? kroui(x[2], p): kronecker(x, p);
     104       13755 :   if (t >= 0) return 0;
     105       21797 :   for (i = lg(L)-1; i; i--)
     106             :   {
     107       14302 :     GEN t = Fp_pow(x, gel(L,i), p);
     108       14302 :     if (equalii(t, p_1) || equali1(t)) return 0;
     109             :   }
     110        7495 :   return 1;
     111             : }
     112             : 
     113             : /* assume p prime, return a generator of all L[i]-Sylows in F_p^*. */
     114             : GEN
     115      357761 : pgener_Fp_local(GEN p, GEN L0)
     116             : {
     117      357761 :   pari_sp av0 = avma;
     118             :   GEN x, p_1, L;
     119      357761 :   if (lgefint(p) == 3)
     120             :   {
     121             :     ulong z;
     122      352872 :     if (p[2] == 2) return gen_1;
     123      257990 :     if (L0) L0 = ZV_to_nv(L0);
     124      257985 :     z = pgener_Fl_local(uel(p,2), L0);
     125      258013 :     return gc_utoipos(av0, z);
     126             :   }
     127        4889 :   p_1 = subiu(p,1); L = is_gener_expo(p, L0);
     128        4890 :   x = utoipos(2);
     129        9899 :   for (;; x[2]++) { if (is_gener_Fp(x, p, p_1, L)) break; }
     130        4890 :   return gc_utoipos(av0, uel(x,2));
     131             : }
     132             : 
     133             : GEN
     134       44134 : pgener_Fp(GEN p) { return pgener_Fp_local(p, NULL); }
     135             : 
     136             : ulong
     137      203274 : pgener_Zl(ulong p)
     138             : {
     139      203274 :   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      203274 :   if (p == 40487) return 10;
     142             : #ifndef LONG_IS_64BIT
     143       29460 :   return pgener_Fl(p);
     144             : #else
     145      173814 :   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      168539 : pgener_Zp(GEN p)
     162             : {
     163      168539 :   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      214579 : rootsof1_Fl(ulong n, ulong p)
     251             : {
     252      214579 :   pari_sp av = avma;
     253      214579 :   GEN L = u_odd_prime_divisors(n); /* 2 implicit in pgener_Fl_local */
     254      214579 :   ulong z = pgener_Fl_local(p, L);
     255      214579 :   z = Fl_powu(z, (p-1) / n, p); /* prim. n-th root of 1 */
     256      214579 :   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   312100188 : ome(long t)
     364             : {
     365   312100188 :   switch(t & 7)
     366             :   {
     367   176905919 :     case 3:
     368   176905919 :     case 5: return 1;
     369   135194269 :     default: return 0;
     370             :   }
     371             : }
     372             : /* t a t_INT, is t = 3,5 mod 8 ? */
     373             : static int
     374     5492364 : gome(GEN t)
     375     5492364 : { return signe(t)? ome( mod2BIL(t) ): 0; }
     376             : 
     377             : /* assume y odd, return kronecker(x,y) * s */
     378             : static long
     379   220107216 : krouu_s(ulong x, ulong y, long s)
     380             : {
     381   220107216 :   ulong x1 = x, y1 = y, z;
     382   999401012 :   while (x1)
     383             :   {
     384   779329853 :     long r = vals(x1);
     385   779401933 :     if (r)
     386             :     {
     387   414825280 :       if (odd(r) && ome(y1)) s = -s;
     388   414717143 :       x1 >>= r;
     389             :     }
     390   779293796 :     if (x1 & y1 & 2) s = -s;
     391   779293796 :     z = y1 % x1; y1 = x1; x1 = z;
     392             :   }
     393   220071159 :   return (y1 == 1)? s: 0;
     394             : }
     395             : 
     396             : long
     397    11612034 : kronecker(GEN x, GEN y)
     398             : {
     399    11612034 :   pari_sp av = avma;
     400    11612034 :   long s = 1, r;
     401             :   ulong xu;
     402             : 
     403    11612034 :   if (typ(x) != t_INT) pari_err_TYPE("kronecker",x);
     404    11612034 :   if (typ(y) != t_INT) pari_err_TYPE("kronecker",y);
     405    11612034 :   switch (signe(y))
     406             :   {
     407          63 :     case -1: y = negi(y); if (signe(x) < 0) s = -1; break;
     408           0 :     case 0: return is_pm1(x);
     409             :   }
     410    11612034 :   r = vali(y);
     411    11612034 :   if (r)
     412             :   {
     413     1347351 :     if (!mpodd(x)) return gc_long(av,0);
     414      322649 :     if (odd(r) && gome(x)) s = -s;
     415      322649 :     y = shifti(y,-r);
     416             :   }
     417    10587331 :   x = modii(x,y);
     418    12969113 :   while (lgefint(x) > 3) /* x < y */
     419             :   {
     420             :     GEN z;
     421     2381831 :     r = vali(x);
     422     2381844 :     if (r)
     423             :     {
     424     1301516 :       if (odd(r) && gome(y)) s = -s;
     425     1301465 :       x = shifti(x,-r);
     426             :     }
     427             :     /* x=3 mod 4 && y=3 mod 4 ? (both are odd here) */
     428     2380393 :     if (mod2BIL(x) & mod2BIL(y) & 2) s = -s;
     429     2380648 :     z = remii(y,x); y = x; x = z;
     430     2381825 :     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    10587282 :   xu = itou(x);
     437    10587286 :   if (!xu) return is_pm1(y)? s: 0;
     438    10521333 :   r = vals(xu);
     439    10521339 :   if (r)
     440             :   {
     441     5600235 :     if (odd(r) && gome(y)) s = -s;
     442     5600235 :     xu >>= r;
     443             :   }
     444             :   /* x=3 mod 4 && y=3 mod 4 ? (both are odd here) */
     445    10521339 :   if (xu & mod2BIL(y) & 2) s = -s;
     446    10521346 :   return gc_long(av, krouu_s(umodiu(y,xu), xu, s));
     447             : }
     448             : 
     449             : long
     450       35973 : krois(GEN x, long y)
     451             : {
     452             :   ulong yu;
     453       35973 :   long s = 1;
     454             : 
     455       35973 :   if (y <= 0)
     456             :   {
     457           7 :     if (y == 0) return is_pm1(x);
     458           0 :     yu = (ulong)-y; if (signe(x) < 0) s = -1;
     459             :   }
     460             :   else
     461       35966 :     yu = (ulong)y;
     462       35966 :   if (!odd(yu))
     463             :   {
     464             :     long r;
     465       14938 :     if (!mpodd(x)) return 0;
     466       11088 :     r = vals(yu); yu >>= r;
     467       11088 :     if (odd(r) && gome(x)) s = -s;
     468             :   }
     469       32116 :   return krouu_s(umodiu(x, yu), yu, s);
     470             : }
     471             : /* assume y != 0 */
     472             : long
     473    27375932 : kroiu(GEN x, ulong y)
     474             : {
     475             :   long r;
     476    27375932 :   if (odd(y)) return krouu_s(umodiu(x,y), y, 1);
     477      299211 :   if (!mpodd(x)) return 0;
     478      204606 :   r = vals(y); y >>= r;
     479      204614 :   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      178605 : krouodd(ulong x, GEN y, long s)
     485             : {
     486             :   long r;
     487      178605 :   if (lgefint(y) == 3) return krouu_s(x, y[2], s);
     488       28776 :   if (!x) return 0; /* y != 1 */
     489       28776 :   r = vals(x);
     490       28776 :   if (r)
     491             :   {
     492       15679 :     if (odd(r) && gome(y)) s = -s;
     493       15679 :     x >>= r;
     494             :   }
     495             :   /* x=3 mod 4 && y=3 mod 4 ? (both are odd here) */
     496       28776 :   if (x & mod2BIL(y) & 2) s = -s;
     497       28776 :   return krouu_s(umodiu(y,x), x, s);
     498             : }
     499             : 
     500             : long
     501      162697 : krosi(long x, GEN y)
     502             : {
     503      162697 :   const pari_sp av = avma;
     504      162697 :   long s = 1, r;
     505      162697 :   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      162697 :   r = vali(y);
     511      162697 :   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      162697 :   if (x < 0) { x = -x; if (mod4(y) == 3) s = -s; }
     518      162697 :   return gc_long(av, krouodd((ulong)x, y, s));
     519             : }
     520             : 
     521             : long
     522       15908 : kroui(ulong x, GEN y)
     523             : {
     524       15908 :   const pari_sp av = avma;
     525       15908 :   long s = 1, r;
     526       15908 :   switch (signe(y))
     527             :   {
     528           0 :     case -1: y = negi(y); break;
     529           0 :     case 0: return x==1UL;
     530             :   }
     531       15908 :   r = vali(y);
     532       15908 :   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       15908 :   return gc_long(av, krouodd(x, y, s));
     539             : }
     540             : 
     541             : long
     542    95816617 : kross(long x, long y)
     543             : {
     544             :   ulong yu;
     545    95816617 :   long s = 1;
     546             : 
     547    95816617 :   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    95747674 :     yu = (ulong)y;
     554    95816589 :   if (!odd(yu))
     555             :   {
     556             :     long r;
     557    23352911 :     if (!odd(x)) return 0;
     558    16438914 :     r = vals(yu); yu >>= r;
     559    16438914 :     if (odd(r) && ome(x)) s = -s;
     560             :   }
     561    88902592 :   x %= (long)yu; if (x < 0) x += yu;
     562    88902592 :   return krouu_s((ulong)x, yu, s);
     563             : }
     564             : 
     565             : long
     566    93217874 : krouu(ulong x, ulong y)
     567             : {
     568             :   long r;
     569    93217874 :   if (odd(y)) return krouu_s(x, y, 1);
     570        5520 :   if (!odd(x)) return 0;
     571        6217 :   r = vals(y); y >>= r;
     572        6217 :   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        7343 : mphilbertoo(GEN x, GEN y)
     581             : {
     582        7343 :   long sx = signe(x), sy = signe(y);
     583        7343 :   if (!sx || !sy) return 0;
     584        7343 :   return (sx < 0 && sy < 0)? -1: 1;
     585             : }
     586             : 
     587             : long
     588      135628 : hilbertii(GEN x, GEN y, GEN p)
     589             : {
     590             :   pari_sp av;
     591             :   long oddvx, oddvy, z;
     592             : 
     593      135628 :   if (!p) return mphilbertoo(x,y);
     594      128306 :   if (is_pm1(p) || signe(p) < 0) pari_err_PRIME("hilbertii",p);
     595      128306 :   if (!signe(x) || !signe(y)) return 0;
     596      128285 :   av = avma;
     597      128285 :   oddvx = odd(Z_pvalrem(x,p,&x));
     598      128285 :   oddvy = odd(Z_pvalrem(y,p,&y));
     599             :   /* x, y are p-units, compute hilbert(x * p^oddvx, y * p^oddvy, p) */
     600      128285 :   if (absequaliu(p, 2))
     601             :   {
     602       12117 :     z = (Mod4(x) == 3 && Mod4(y) == 3)? -1: 1;
     603       12117 :     if (oddvx && gome(y)) z = -z;
     604       12117 :     if (oddvy && gome(x)) z = -z;
     605             :   }
     606             :   else
     607             :   {
     608      116168 :     z = (oddvx && oddvy && mod4(p) == 3)? -1: 1;
     609      116168 :     if (oddvx && kronecker(y,p) < 0) z = -z;
     610      116168 :     if (oddvy && kronecker(x,p) < 0) z = -z;
     611             :   }
     612      128285 :   return gc_long(av, z);
     613             : }
     614             : 
     615             : static void
     616         196 : err_prec(void) { pari_err_PREC("hilbert"); }
     617             : static void
     618         161 : err_p(GEN p, GEN q) { pari_err_MODULUS("hilbert", p,q); }
     619             : static void
     620          56 : err_oo(GEN p) { pari_err_MODULUS("hilbert", p, strtoGENstr("oo")); }
     621             : 
     622             : /* x t_INTMOD, *pp = prime or NULL [ unset, set it to x.mod ].
     623             :  * Return lift(x) provided it's p-adic accuracy is large enough to decide
     624             :  * hilbert()'s value [ problem at p = 2 ] */
     625             : static GEN
     626         420 : lift_intmod(GEN x, GEN *pp)
     627             : {
     628         420 :   GEN p = *pp, N = gel(x,1);
     629         420 :   x = gel(x,2);
     630         420 :   if (!p)
     631             :   {
     632         266 :     *pp = p = N;
     633         266 :     switch(itos_or_0(p))
     634             :     {
     635         126 :       case 2:
     636         126 :       case 4: err_prec();
     637             :     }
     638         140 :     return x;
     639             :   }
     640         154 :   if (!signe(p)) err_oo(N);
     641         112 :   if (absequaliu(p,2))
     642          42 :   { if (vali(N) <= 2) err_prec(); }
     643             :   else
     644          70 :   { if (!dvdii(N,p)) err_p(N,p); }
     645          28 :   if (!signe(x)) err_prec();
     646          21 :   return x;
     647             : }
     648             : /* x t_PADIC, *pp = prime or NULL [ unset, set it to x.p ].
     649             :  * Return lift(x)*p^(v(x) mod 2) provided it's p-adic accuracy is large enough
     650             :  * to decide hilbert()'s value [ problem at p = 2 ]*/
     651             : static GEN
     652         210 : lift_padic(GEN x, GEN *pp)
     653             : {
     654         210 :   GEN p = *pp, q = gel(x,2), y = gel(x,4);
     655         210 :   if (!p) *pp = p = q;
     656         147 :   else if (!equalii(p,q)) err_p(p, q);
     657         105 :   if (absequaliu(p,2) && precp(x) <= 2) err_prec();
     658          70 :   if (!signe(y)) err_prec();
     659          70 :   return odd(valp(x))? mulii(p,y): y;
     660             : }
     661             : 
     662             : long
     663       60665 : hilbert(GEN x, GEN y, GEN p)
     664             : {
     665       60665 :   pari_sp av = avma;
     666       60665 :   long tx = typ(x), ty = typ(y);
     667             : 
     668       60665 :   if (p && typ(p) != t_INT) pari_err_TYPE("hilbert",p);
     669       60665 :   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       60588 :   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       60574 :   if (tx == t_INTMOD) { x = lift_intmod(x, &p); tx = t_INT; }
     692       60371 :   if (ty == t_INTMOD) { y = lift_intmod(y, &p); ty = t_INT; }
     693             : 
     694       60315 :   if (tx == t_PADIC) { x = lift_padic(x, &p); tx = t_INT; }
     695       60252 :   if (ty == t_PADIC) { y = lift_padic(y, &p); ty = t_INT; }
     696             : 
     697       60175 :   if (tx == t_FRAC) { tx = t_INT; x = p? mulii(gel(x,1),gel(x,2)): gel(x,1); }
     698       60175 :   if (ty == t_FRAC) { ty = t_INT; y = p? mulii(gel(y,1),gel(y,2)): gel(y,1); }
     699             : 
     700       60175 :   if (tx != t_INT || ty != t_INT) pari_err_TYPE2("hilbert",x,y);
     701       60175 :   if (p && !signe(p)) p = NULL;
     702       60175 :   return gc_long(av, hilbertii(x,y,p));
     703             : }
     704             : 
     705             : /*******************************************************************/
     706             : /*                       SQUARE ROOT MODULO p                      */
     707             : /*******************************************************************/
     708             : static void
     709     2078507 : checkp(ulong q, ulong p)
     710     2078507 : { 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    11078356 : nonsquare1_Fl(ulong p)
     714             : {
     715             :   forprime_t S;
     716             :   ulong q;
     717    11078356 :   if ((p & 7UL) != 1) return 2UL;
     718     4104123 :   q = p % 3; if (q == 2) return 3UL;
     719     1303148 :   checkp(q, p);
     720     1309710 :   q = p % 5; if (q == 2 || q == 3) return 5UL;
     721      482051 :   checkp(q, p);
     722      482051 :   q = p % 7; if (q != 4 && q >= 3) return 7UL;
     723      176525 :   checkp(q, p);
     724             :   /* log^2(2^64) < 1968 is enough under GRH (and p^(1/4)log(p) without it)*/
     725      176547 :   u_forprime_init(&S, 11, 1967);
     726      286988 :   while ((q = u_forprime_next(&S)))
     727             :   {
     728      286986 :     if (krouu(q, p) < 0) return q;
     729      110435 :     checkp(q, p);
     730             :   }
     731           0 :   checkp(0, p);
     732             :   return 0; /*LCOV_EXCL_LINE*/
     733             : }
     734             : /* p > 2 a prime */
     735             : ulong
     736        7893 : nonsquare_Fl(ulong p)
     737        7893 : { return ((p & 3UL) == 3)? p-1: nonsquare1_Fl(p); }
     738             : 
     739             : /* allow pi = 0 */
     740             : ulong
     741      173626 : Fl_2gener_pre(ulong p, ulong pi)
     742             : {
     743      173626 :   ulong p1 = p-1;
     744      173626 :   long e = vals(p1);
     745      173614 :   if (e == 1) return p1;
     746       64698 :   return Fl_powu_pre(nonsquare1_Fl(p), p1 >> e, p, pi);
     747             : }
     748             : 
     749             : ulong
     750       67793 : Fl_2gener_pre_i(ulong  ns, ulong p, ulong pi)
     751             : {
     752       67793 :   ulong p1 = p-1;
     753       67793 :   long e = vals(p1);
     754       67793 :   if (e == 1) return p1;
     755       26020 :   return Fl_powu_pre(ns, p1 >> e, p, pi);
     756             : }
     757             : 
     758             : static ulong
     759    11495909 : Fl_sqrt_i(ulong a, ulong y, ulong p)
     760             : {
     761             :   long i, e, k;
     762             :   ulong p1, q, v, w;
     763             : 
     764    11495909 :   if (!a) return 0;
     765    10243682 :   p1 = p - 1; e = vals(p1);
     766    10243450 :   if (e == 0) /* p = 2 */
     767             :   {
     768      421559 :     if (p != 2) pari_err_PRIME("Fl_sqrt [modulus]",utoi(p));
     769      422606 :     return ((a & 1) == 0)? 0: 1;
     770             :   }
     771     9821891 :   if (e == 1)
     772             :   {
     773     4728353 :     v = Fl_powu(a, (p+1) >> 2, p);
     774     4728542 :     if (Fl_sqr(v, p) != a) return ~0UL;
     775     4723699 :     p1 = p - v; if (v > p1) v = p1;
     776     4723699 :     return v;
     777             :   }
     778     5093538 :   q = p1 >> e; /* q = (p-1)/2^oo is odd */
     779     5093538 :   p1 = Fl_powu(a, q >> 1, p); /* a ^ [(q-1)/2] */
     780     5093632 :   if (!p1) return 0;
     781     5093632 :   v = Fl_mul(a, p1, p);
     782     5093646 :   w = Fl_mul(v, p1, p);
     783     5093738 :   if (!y) y = Fl_powu(nonsquare1_Fl(p), q, p);
     784     8663633 :   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     3571826 :     p1 = Fl_sqr(w, p);
     788     5881970 :     for (k=1; p1 != 1 && k < e; k++) p1 = Fl_sqr(p1, p);
     789     3571833 :     if (k == e) return ~0UL;
     790             :     /* w ^ (2^k) = 1 --> w = y ^ (u * 2^(e-k)), u odd */
     791     3569756 :     p1 = y;
     792     4688866 :     for (i=1; i < e-k; i++) p1 = Fl_sqr(p1, p);
     793     3569756 :     y = Fl_sqr(p1, p); e = k;
     794     3569781 :     w = Fl_mul(y, w, p);
     795     3569815 :     v = Fl_mul(v, p1, p);
     796             :   }
     797     5091807 :   p1 = p - v; if (v > p1) v = p1;
     798     5091807 :   return v;
     799             : }
     800             : 
     801             : /* Tonelli-Shanks. Assume p is prime and (a,p) != -1. Allow pi = 0 */
     802             : ulong
     803    32290597 : 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    32290597 :   if (!pi) return Fl_sqrt_i(a, y, p);
     809    20794801 :   if (!a) return 0;
     810    20670902 :   p1 = p - 1; e = vals(p1);
     811    20678218 :   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    20689568 :   if (e == 1)
     817             :   {
     818    14716161 :     v = Fl_powu_pre(a, (p+1) >> 2, p, pi);
     819    14672229 :     if (Fl_sqr_pre(v, p, pi) != a) return ~0UL;
     820    14686448 :     p1 = p - v; if (v > p1) v = p1;
     821    14686448 :     return v;
     822             :   }
     823     5973407 :   q = p1 >> e; /* q = (p-1)/2^oo is odd */
     824     5973407 :   p1 = Fl_powu_pre(a, q >> 1, p, pi); /* a ^ [(q-1)/2] */
     825     5968696 :   if (!p1) return 0;
     826     5968696 :   v = Fl_mul_pre(a, p1, p, pi);
     827     5969926 :   w = Fl_mul_pre(v, p1, p, pi);
     828     5969577 :   if (!y) y = Fl_powu_pre(nonsquare1_Fl(p), q, p, pi);
     829    11362986 :   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     5393988 :     p1 = Fl_sqr_pre(w,p,pi);
     833    10137092 :     for (k=1; p1 != 1 && k < e; k++) p1 = Fl_sqr_pre(p1,p,pi);
     834     5395855 :     if (k == e) return ~0UL;
     835             :     /* w ^ (2^k) = 1 --> w = y ^ (u * 2^(e-k)), u odd */
     836     5395763 :     p1 = y;
     837     7068512 :     for (i=1; i < e-k; i++) p1 = Fl_sqr_pre(p1, p, pi);
     838     5395769 :     y = Fl_sqr_pre(p1, p, pi); e = k;
     839     5396146 :     w = Fl_mul_pre(y, w, p, pi);
     840     5395690 :     v = Fl_mul_pre(v, p1, p, pi);
     841             :   }
     842     5968998 :   p1 = p - v; if (v > p1) v = p1;
     843     5968998 :   return v;
     844             : }
     845             : 
     846             : ulong
     847    11511734 : Fl_sqrt(ulong a, ulong p)
     848    11511734 : { ulong pi = (p & HIGHMASK)? get_Fl_red(p): 0; return Fl_sqrt_pre_i(a, 0, p, pi); }
     849             : 
     850             : ulong
     851    20622190 : Fl_sqrt_pre(ulong a, ulong p, ulong pi)
     852    20622190 : { return Fl_sqrt_pre_i(a, 0, p, pi); }
     853             : 
     854             : /* allow pi = 0 */
     855             : static ulong
     856       51035 : Fl_lgener_pre_all(ulong l, long e, ulong r, ulong p, ulong pi, ulong *pt_m)
     857             : {
     858       51035 :   ulong x, y, m, le1 = upowuu(l, e-1);
     859       51034 :   for (x = 2; ; x++)
     860             :   {
     861       79000 :     y = Fl_powu_pre(x, r, p, pi);
     862       78998 :     if (y==1) continue;
     863       62341 :     m = Fl_powu_pre(y, le1, p, pi);
     864       62346 :     if (m != 1) break;
     865             :   }
     866       51037 :   *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      124189 : 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      124189 :   if (a==0) return a;
     879      124183 :   u2 = Fl_inv(l%r, r);
     880      124184 :   v = Fl_powu_pre(a, u2, p, pi);
     881      124182 :   w = Fl_powu_pre(v, l, p, pi);
     882      124177 :   w = pi? Fl_mul_pre(w, Fl_inv(a, p), p, pi): Fl_div(w, a, p);
     883      124166 :   if (w==1) return v;
     884       49684 :   if (y==0) y = Fl_lgener_pre_all(l, e, r, p, pi, &m);
     885       72279 :   while (w!=1)
     886             :   {
     887       54923 :     ulong k = 0;
     888       54923 :     p1 = w;
     889             :     do
     890             :     {
     891       81422 :       z = p1; p1 = Fl_powu_pre(p1, l, p, pi);
     892       81422 :       if (++k == e) return ULONG_MAX;
     893       49092 :     } while (p1!=1);
     894       22593 :     dl = Fl_log_pre(z, m, l, p, pi);
     895       22593 :     dl = Fl_neg(dl, l);
     896       22593 :     p1 = Fl_powu_pre(y,dl*upowuu(l,e-k-1),p,pi);
     897       22593 :     m = Fl_powu_pre(m, dl, p, pi);
     898       22593 :     e = k;
     899       22593 :     v = pi? Fl_mul_pre(p1,v,p,pi): Fl_mul(p1,v,p);
     900       22593 :     y = Fl_powu_pre(p1,l,p,pi);
     901       22593 :     w = pi? Fl_mul_pre(y,w,p,pi): Fl_mul(y,w,p);
     902             :   }
     903       17356 :   return v;
     904             : }
     905             : 
     906             : /* allow pi = 0 */
     907             : static ulong
     908      121855 : Fl_sqrtl_i(ulong a, ulong l, ulong p, ulong pi, ulong y, ulong m)
     909             : {
     910      121855 :   ulong r, e = u_lvalrem(p-1, l, &r);
     911      121856 :   return Fl_sqrtl_raw(a, l, e, r, p, pi, y, m);
     912             : }
     913             : /* allow pi = 0 */
     914             : ulong
     915      121854 : Fl_sqrtl_pre(ulong a, ulong l, ulong p, ulong pi)
     916      121854 : { 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       85685 : Fl_sqrtn_pre(ulong a, long n, ulong p, ulong pi, ulong *zetan)
     926             : {
     927       85685 :   ulong m, q = p-1, z;
     928       85685 :   ulong nn = n >= 0 ? (ulong)n: -(ulong)n;
     929       85685 :   if (a==0)
     930             :   {
     931       48153 :     if (n < 0) pari_err_INV("Fl_sqrtn", mkintmod(gen_0,utoi(p)));
     932       48146 :     if (zetan) *zetan = 1UL;
     933       48146 :     return 0;
     934             :   }
     935       37532 :   if (n==1)
     936             :   {
     937         420 :     if (zetan) *zetan = 1;
     938         420 :     return n < 0? Fl_inv(a,p): a;
     939             :   }
     940       37112 :   if (n==2)
     941             :   {
     942        5710 :     if (zetan) *zetan = p-1;
     943        5710 :     return Fl_sqrt_pre_i(a, 0, p, pi);
     944             :   }
     945       31402 :   if (a == 1 && !zetan) return a;
     946        7973 :   m = ugcd(nn, q);
     947        7973 :   z = 1;
     948        7973 :   if (m!=1)
     949             :   {
     950        1400 :     GEN F = factoru(m);
     951             :     long i, j, e;
     952             :     ulong r, zeta, y, l;
     953        3108 :     for (i = nbrows(F); i; i--)
     954             :     {
     955        1771 :       l = ucoeff(F,i,1);
     956        1771 :       j = ucoeff(F,i,2);
     957        1771 :       e = u_lvalrem(q,l, &r);
     958        1771 :       y = Fl_lgener_pre_all(l, e, r, p, pi, &zeta);
     959        1771 :       if (zetan)
     960             :       {
     961         112 :         ulong Y = Fl_powu_pre(y, upowuu(l,e-j), p, pi);
     962         112 :         z = pi? Fl_mul_pre(z, Y, p, pi): Fl_mul(z, Y, p);
     963             :       }
     964        1771 :       if (a!=1)
     965             :         do
     966             :         {
     967        2331 :           a = Fl_sqrtl_raw(a, l, e, r, p, pi, y, zeta);
     968        2317 :           if (a==ULONG_MAX) return ULONG_MAX;
     969        2268 :         } while (--j);
     970             :     }
     971             :   }
     972        7910 :   if (m != nn)
     973             :   {
     974        6594 :     ulong qm = q/m, nm = (nn/m) % qm;
     975        6594 :     a = Fl_powu_pre(a, Fl_inv(nm, qm), p, pi);
     976             :   }
     977        7910 :   if (n < 0) a = Fl_inv(a, p);
     978        7910 :   if (zetan) *zetan = z;
     979        7910 :   return a;
     980             : }
     981             : 
     982             : ulong
     983       85685 : Fl_sqrtn(ulong a, long n, ulong p, ulong *zetan)
     984             : {
     985       85685 :   ulong pi = (p & HIGHMASK)? get_Fl_red(p): 0;
     986       85685 :   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         449 : sqrt_Cipolla_sqr(void *data, GEN y)
    1002             : {
    1003         449 :   GEN u = gel(y,1), v = gel(y,2), p = gel(data,2), n = gel(data,3);
    1004         449 :   GEN u2 = sqri(u), v2 = sqri(v);
    1005         449 :   v = subii(sqri(addii(v,u)), addii(u2,v2));
    1006         449 :   u = addii(u2, mulii(v2,n));
    1007         449 :   retmkvec2(modii(u,p), modii(v,p));
    1008             : }
    1009             : /* compute (t+X) y^2 */
    1010             : static GEN
    1011          23 : sqrt_Cipolla_msqr(void *data, GEN y)
    1012             : {
    1013          23 :   GEN u = gel(y,1), v = gel(y,2), a = gel(data,1), p = gel(data,2);
    1014          23 :   ulong t = gel(data,4)[2];
    1015          23 :   GEN d = addii(u, mului(t,v)), d2 = sqri(d);
    1016          23 :   GEN b = remii(mulii(a,v), p);
    1017          23 :   u = subii(mului(t,d2), mulii(b,addii(u,d)));
    1018          23 :   v = subii(d2, mulii(b,v));
    1019          23 :   retmkvec2(modii(u,p), modii(v,p));
    1020             : }
    1021             : /* assume a reduced mod p [ otherwise correct but inefficient ] */
    1022             : static GEN
    1023           8 : sqrt_Cipolla(GEN a, GEN p)
    1024             : {
    1025             :   pari_sp av;
    1026             :   GEN u, v, n, y, pov2;
    1027             :   ulong t;
    1028             : 
    1029           8 :   if (kronecker(a, p) < 0) return NULL;
    1030           8 :   pov2 = shifti(p,-1); /* center to avoid multiplying by huge base*/
    1031           8 :   if (cmpii(a,pov2) > 0) a = subii(a,p);
    1032           8 :   av = avma;
    1033          41 :   for (t=1; ; t++, set_avma(av))
    1034             :   {
    1035          41 :     n = subsi((long)(t*t), a);
    1036          41 :     if (kronecker(n, p) < 0) break;
    1037             :   }
    1038             : 
    1039             :   /* compute (t+X)^((p-1)/2) =: u+vX */
    1040           8 :   u = utoipos(t);
    1041           8 :   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             : 
    1050           8 :   v = Fp_mul(gel(y, 2), a, p);
    1051           8 :   if (cmpii(v,pov2) > 0) v = subii(p,v);
    1052           8 :   return v;
    1053             : }
    1054             : 
    1055             : /* Return NULL if p is found to be composite */
    1056             : static GEN
    1057        6097 : Fp_2gener_all(long e, GEN p)
    1058             : {
    1059             :   GEN y, m;
    1060             :   long k;
    1061        6097 :   GEN q = shifti(subiu(p,1), -e); /* q = (p-1)/2^oo is odd */
    1062        6097 :   if (e==0 && !equaliu(p,2)) return NULL;
    1063        6097 :   for (k=2; ; k++)
    1064       12336 :   {
    1065       18433 :     long i = krosi(k, p);
    1066       18433 :     if (i >= 0)
    1067             :     {
    1068       12336 :       if (i) continue;
    1069           0 :       return NULL;
    1070             :     }
    1071        6097 :     y = m = Fp_pow(utoi(k), q, p);
    1072       17986 :     for (i=1; i<e; i++)
    1073       11889 :       if (equali1(m = Fp_sqr(m, p))) break;
    1074        6097 :     if (i == e) break; /* success */
    1075             :   }
    1076        6097 :   return y;
    1077             : }
    1078             : 
    1079             : /* Return NULL if p is found to be composite */
    1080             : GEN
    1081        3185 : Fp_2gener(GEN p)
    1082        3185 : { return Fp_2gener_all(vali(subis(p,1)),p); }
    1083             : 
    1084             : GEN
    1085       19931 : Fp_2gener_i(GEN ns, GEN p)
    1086             : {
    1087       19931 :   GEN p1 = subiu(p,1);
    1088       19931 :   long e = vali(p1);
    1089       19931 :   if (e == 1) return p1;
    1090       18546 :   return Fp_pow(ns, shifti(p1,-e), p);
    1091             : }
    1092             : 
    1093             : /* smallest square root */
    1094             : static GEN
    1095       78327 : choose_sqrt(GEN v, GEN p)
    1096             : {
    1097       78327 :   pari_sp av = avma;
    1098       78327 :   GEN q = subii(p,v);
    1099       78323 :   if (cmpii(v,q) > 0) v = q; else set_avma(av);
    1100       78325 :   return v;
    1101             : }
    1102             : /* Tonelli-Shanks. Assume p is prime and return NULL if (a,p) = -1. */
    1103             : GEN
    1104     3382295 : Fp_sqrt_i(GEN a, GEN y, GEN p)
    1105             : {
    1106     3382295 :   pari_sp av = avma;
    1107             :   long i, k, e;
    1108             :   GEN p1, q, v, w;
    1109             : 
    1110     3382295 :   if (lgefint(p) == 3)
    1111             :   {
    1112     3303454 :     ulong pp = uel(p,2), u = umodiu(a, pp);
    1113     3303487 :     if (!u) return gen_0;
    1114     3006980 :     u = Fl_sqrt(u, pp);
    1115     3007022 :     return (u == ~0UL)? NULL: utoipos(u);
    1116             :   }
    1117             : 
    1118       78841 :   a = modii(a, p); if (!signe(a)) return gc_const(av, gen_0);
    1119       78720 :   p1 = subiu(p,1); e = vali(p1);
    1120       78738 :   if (e <= 2)
    1121             :   { /* direct formulas more efficient */
    1122             :     pari_sp av2;
    1123       38392 :     if (e == 0) pari_err_PRIME("Fp_sqrt [modulus]",p); /* p != 2 */
    1124       38394 :     if (e == 1)
    1125             :     {
    1126       19638 :       q = addiu(shifti(p1,-2),1); /* (p+1) / 4 */
    1127       19627 :       v = Fp_pow(a, q, p);
    1128             :     }
    1129             :     else
    1130             :     { /* Atkin's formula */
    1131       18756 :       GEN i, a2 = shifti(a,1);
    1132       18750 :       if (cmpii(a2,p) >= 0) a2 = subii(a2,p);
    1133       18755 :       q = shifti(p1, -3); /* (p-5)/8 */
    1134       18757 :       v = Fp_pow(a2, q, p);
    1135       18760 :       i = Fp_mul(a2, Fp_sqr(v,p), p); /* i^2 = -1 */
    1136       18760 :       v = Fp_mul(a, Fp_mul(v, subiu(i,1), p), p);
    1137             :     }
    1138       38418 :     av2 = avma;
    1139             :     /* must check equality in case (a/p) = -1 or p not prime */
    1140       38418 :     e = equalii(Fp_sqr(v,p), a); set_avma(av2);
    1141       38418 :     return e? gerepileuptoint(av,choose_sqrt(v,p)): NULL;
    1142             :   }
    1143             :   /* On average, Cipolla is better than Tonelli/Shanks if and only if
    1144             :    * e(e-1) > 8*log2(n)+20, see LNCS 2286 pp 430 [GTL] */
    1145       40346 :   if (e*(e-1) > 20 + 8 * expi(p))
    1146             :   {
    1147           8 :     v = sqrt_Cipolla(a,p); if (!v) return gc_NULL(av);
    1148           8 :     return gerepileuptoint(av,v);
    1149             :   }
    1150       40339 :   if (!y)
    1151             :   {
    1152        2912 :     y = Fp_2gener_all(e, p);
    1153        2912 :     if (!y) pari_err_PRIME("Fp_sqrt [modulus]",p);
    1154             :   }
    1155       40339 :   q = shifti(p1,-e); /* q = (p-1)/2^oo is odd */
    1156       40338 :   p1 = Fp_pow(a, shifti(q,-1), p); /* a ^ (q-1)/2 */
    1157       40342 :   v = Fp_mul(a, p1, p);
    1158       40343 :   w = Fp_mul(v, p1, p);
    1159       96095 :   while (!equali1(w))
    1160             :   { /* a*w = v^2, y primitive 2^e-th root of 1
    1161             :        a square --> w even power of y, hence w^(2^(e-1)) = 1 */
    1162       55796 :     p1 = Fp_sqr(w,p);
    1163      115050 :     for (k=1; !equali1(p1) && k < e; k++) p1 = Fp_sqr(p1,p);
    1164       55791 :     if (k == e) return gc_NULL(av); /* p composite or (a/p) != 1 */
    1165             :     /* w ^ (2^k) = 1 --> w = y ^ (u * 2^(e-k)), u odd */
    1166       55748 :     p1 = y;
    1167       79619 :     for (i=1; i < e-k; i++) p1 = Fp_sqr(p1,p);
    1168       55748 :     y = Fp_sqr(p1, p); e = k;
    1169       55754 :     w = Fp_mul(y, w, p);
    1170       55752 :     v = Fp_mul(v, p1, p);
    1171       55753 :     if (gc_needed(av,1))
    1172             :     {
    1173           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"Fp_sqrt");
    1174           0 :       gerepileall(av,3, &y,&w,&v);
    1175             :     }
    1176             :   }
    1177       40298 :   return gerepileuptoint(av, choose_sqrt(v,p));
    1178             : }
    1179             : 
    1180             : GEN
    1181     3323902 : Fp_sqrt(GEN a, GEN p)
    1182     3323902 : { return Fp_sqrt_i(a, NULL, p); }
    1183             : 
    1184             : /*********************************************************************/
    1185             : /**                        GCD & BEZOUT                             **/
    1186             : /*********************************************************************/
    1187             : 
    1188             : GEN
    1189    43641036 : lcmii(GEN x, GEN y)
    1190             : {
    1191             :   pari_sp av;
    1192             :   GEN a, b;
    1193    43641036 :   if (!signe(x) || !signe(y)) return gen_0;
    1194    43641116 :   av = avma; a = gcdii(x,y);
    1195    43639449 :   if (absequalii(a,y)) { set_avma(av); return absi(x); }
    1196     8639931 :   if (!equali1(a)) y = diviiexact(y,a);
    1197     8639946 :   b = mulii(x,y); setabssign(b); return gerepileuptoint(av, b);
    1198             : }
    1199             : 
    1200             : /* given x in assume 0 < x < N; return u in (Z/NZ)^* such that u x = gcd(x,N) (mod N);
    1201             :  * set *pd = gcd(x,N) */
    1202             : GEN
    1203     6610419 : Fp_invgen(GEN x, GEN N, GEN *pd)
    1204             : {
    1205             :   GEN d, d0, e, v;
    1206     6610419 :   if (lgefint(N) == 3)
    1207             :   {
    1208     5829200 :     ulong dd, NN = N[2], xx = umodiu(x,NN);
    1209     5829760 :     if (!xx) { *pd = N; return gen_0; }
    1210     5829760 :     xx = Fl_invgen(xx, NN, &dd);
    1211     5830378 :     *pd = utoi(dd); return utoi(xx);
    1212             :   }
    1213      781219 :   *pd = d = bezout(x, N, &v, NULL);
    1214      781232 :   if (equali1(d)) return v;
    1215             :   /* vx = gcd(x,N) (mod N), v coprime to N/d but need not be coprime to N */
    1216      691263 :   e = diviiexact(N,d);
    1217      691263 :   d0 = Z_ppo(d, e); /* d = d0 d1, d0 coprime to N/d, rad(d1) | N/d */
    1218      691263 :   if (equali1(d0)) return v;
    1219      527976 :   if (!equalii(d,d0)) e = lcmii(e, diviiexact(d,d0));
    1220      527976 :   return Z_chinese_coprime(v, gen_1, e, d0, mulii(e,d0));
    1221             : }
    1222             : 
    1223             : /*********************************************************************/
    1224             : /**                      CHINESE REMAINDERS                         **/
    1225             : /*********************************************************************/
    1226             : 
    1227             : /* Chinese Remainder Theorem.  x and y must have the same type (integermod,
    1228             :  * polymod, or polynomial/vector/matrix recursively constructed with these
    1229             :  * as coefficients). Creates (with the same type) a z in the same residue
    1230             :  * class as x and the same residue class as y, if it is possible.
    1231             :  *
    1232             :  * We also allow (during recursion) two identical objects even if they are
    1233             :  * not integermod or polymod. For example:
    1234             :  *
    1235             :  * ? x = [1, Mod(5, 11), Mod(X + Mod(2, 7), X^2 + 1)];
    1236             :  * ? y = [1, Mod(7, 17), Mod(X + Mod(0, 3), X^2 + 1)];
    1237             :  * ? chinese(x, y)
    1238             :  * %3 = [1, Mod(16, 187), Mod(X + mod(9, 21), X^2 + 1)] */
    1239             : 
    1240             : static GEN
    1241     2401623 : gen_chinese(GEN x, GEN(*f)(GEN,GEN))
    1242             : {
    1243     2401623 :   GEN z = gassoc_proto(f,x,NULL);
    1244     2401616 :   if (z == gen_1) retmkintmod(gen_0,gen_1);
    1245     2401581 :   return z;
    1246             : }
    1247             : 
    1248             : /* x t_INTMOD, y t_POLMOD; promote x to t_POLMOD mod Pol(x.mod) then
    1249             :  * call chinese: makes Mod(0,1) a better "neutral" element */
    1250             : static GEN
    1251          21 : chinese_intpol(GEN x,GEN y)
    1252             : {
    1253          21 :   pari_sp av = avma;
    1254          21 :   GEN z = mkpolmod(gel(x,2), scalarpol_shallow(gel(x,1), varn(gel(y,1))));
    1255          21 :   return gerepileupto(av, chinese(z, y));
    1256             : }
    1257             : 
    1258             : GEN
    1259        2345 : chinese1(GEN x) { return gen_chinese(x,chinese); }
    1260             : 
    1261             : GEN
    1262        5432 : chinese(GEN x, GEN y)
    1263             : {
    1264             :   pari_sp av;
    1265        5432 :   long tx = typ(x), ty;
    1266             :   GEN z,p1,p2,d,u,v;
    1267             : 
    1268        5432 :   if (!y) return chinese1(x);
    1269        5383 :   if (gidentical(x,y)) return gcopy(x);
    1270        5376 :   ty = typ(y);
    1271        5376 :   if (tx == ty) switch(tx)
    1272             :   {
    1273        3829 :     case t_POLMOD:
    1274             :     {
    1275        3829 :       GEN A = gel(x,1), B = gel(y,1);
    1276        3829 :       GEN a = gel(x,2), b = gel(y,2);
    1277        3829 :       if (varn(A)!=varn(B)) pari_err_VAR("chinese",A,B);
    1278        3829 :       if (RgX_equal(A,B)) retmkpolmod(chinese(a,b), gcopy(A)); /*same modulus*/
    1279        3829 :       av = avma;
    1280        3829 :       d = RgX_extgcd(A,B,&u,&v);
    1281        3829 :       p2 = gsub(b, a);
    1282        3829 :       if (!gequal0(gmod(p2, d))) break;
    1283        3829 :       p1 = gdiv(A,d);
    1284        3829 :       p2 = gadd(a, gmul(gmul(u,p1), p2));
    1285             : 
    1286        3829 :       z = cgetg(3, t_POLMOD);
    1287        3829 :       gel(z,1) = gmul(p1,B);
    1288        3829 :       gel(z,2) = gmod(p2,gel(z,1));
    1289        3829 :       return gerepileupto(av, z);
    1290             :     }
    1291        1505 :     case t_INTMOD:
    1292             :     {
    1293        1505 :       GEN A = gel(x,1), B = gel(y,1);
    1294        1505 :       GEN a = gel(x,2), b = gel(y,2), c, d, C, U;
    1295        1505 :       z = cgetg(3,t_INTMOD);
    1296        1505 :       Z_chinese_pre(A, B, &C, &U, &d);
    1297        1505 :       c = Z_chinese_post(a, b, C, U, d);
    1298        1505 :       if (!c) pari_err_OP("chinese", x,y);
    1299        1505 :       set_avma((pari_sp)z);
    1300        1505 :       gel(z,1) = icopy(C);
    1301        1505 :       gel(z,2) = icopy(c); return z;
    1302             :     }
    1303          14 :     case t_POL:
    1304             :     {
    1305          14 :       long i, lx = lg(x), ly = lg(y);
    1306          14 :       if (varn(x) != varn(y)) break;
    1307          14 :       if (lx < ly) { swap(x,y); lswap(lx,ly); }
    1308          14 :       z = cgetg(lx, t_POL); z[1] = x[1];
    1309          42 :       for (i=2; i<ly; i++) gel(z,i) = chinese(gel(x,i),gel(y,i));
    1310          14 :       if (i < lx)
    1311             :       {
    1312          14 :         GEN _0 = Rg_get_0(y);
    1313          28 :         for (   ; i<lx; i++) gel(z,i) = chinese(gel(x,i),_0);
    1314             :       }
    1315          14 :       return z;
    1316             :     }
    1317           7 :     case t_VEC: case t_COL: case t_MAT:
    1318             :     {
    1319             :       long i, lx;
    1320           7 :       z = cgetg_copy(x, &lx); if (lx!=lg(y)) break;
    1321          21 :       for (i=1; i<lx; i++) gel(z,i) = chinese(gel(x,i),gel(y,i));
    1322           7 :       return z;
    1323             :     }
    1324             :   }
    1325          21 :   if (tx == t_POLMOD && ty == t_INTMOD) return chinese_intpol(y,x);
    1326           7 :   if (ty == t_POLMOD && tx == t_INTMOD) return chinese_intpol(x,y);
    1327           0 :   pari_err_OP("chinese",x,y);
    1328             :   return NULL; /* LCOV_EXCL_LINE */
    1329             : }
    1330             : 
    1331             : /* init chinese(Mod(.,A), Mod(.,B)) */
    1332             : void
    1333      267283 : Z_chinese_pre(GEN A, GEN B, GEN *pC, GEN *pU, GEN *pd)
    1334             : {
    1335      267283 :   GEN u, d = bezout(A,B,&u,NULL); /* U = u(A/d), u(A/d) + v(B/d) = 1 */
    1336      267285 :   GEN t = diviiexact(A,d);
    1337      267262 :   *pU = mulii(u, t);
    1338      267265 :   *pC = mulii(t, B);
    1339      267257 :   if (pd) *pd = d;
    1340      267257 : }
    1341             : /* Assume C = lcm(A, B), U = 0 mod (A/d), U = 1 mod (B/d), a = b mod d,
    1342             :  * where d = gcd(A,B) or NULL, return x = a (mod A), b (mod B).
    1343             :  * If d not NULL, check whether a = b mod d. */
    1344             : GEN
    1345     2981548 : Z_chinese_post(GEN a, GEN b, GEN C, GEN U, GEN d)
    1346             : {
    1347             :   GEN b_a;
    1348     2981548 :   if (!signe(a))
    1349             :   {
    1350      788655 :     if (d && !dvdii(b, d)) return NULL;
    1351      788655 :     return Fp_mul(b, U, C);
    1352             :   }
    1353     2192893 :   b_a = subii(b,a);
    1354     2192893 :   if (d && !dvdii(b_a, d)) return NULL;
    1355     2192893 :   return modii(addii(a, mulii(U, b_a)), C);
    1356             : }
    1357             : static ulong
    1358     2178065 : u_chinese_post(ulong a, ulong b, ulong C, ulong U)
    1359             : {
    1360     2178065 :   if (!a) return Fl_mul(b, U, C);
    1361     2178065 :   return Fl_add(a, Fl_mul(U, Fl_sub(b,a,C), C), C);
    1362             : }
    1363             : 
    1364             : GEN
    1365        2142 : Z_chinese(GEN a, GEN b, GEN A, GEN B)
    1366             : {
    1367        2142 :   pari_sp av = avma;
    1368        2142 :   GEN C, U; Z_chinese_pre(A, B, &C, &U, NULL);
    1369        2142 :   return gerepileuptoint(av, Z_chinese_post(a,b, C, U, NULL));
    1370             : }
    1371             : GEN
    1372      263579 : Z_chinese_all(GEN a, GEN b, GEN A, GEN B, GEN *pC)
    1373             : {
    1374      263579 :   GEN U; Z_chinese_pre(A, B, pC, &U, NULL);
    1375      263554 :   return Z_chinese_post(a,b, *pC, U, NULL);
    1376             : }
    1377             : 
    1378             : /* return lift(chinese(a mod A, b mod B))
    1379             :  * assume(A,B)=1, a,b,A,B integers and C = A*B */
    1380             : GEN
    1381      529236 : Z_chinese_coprime(GEN a, GEN b, GEN A, GEN B, GEN C)
    1382             : {
    1383      529236 :   pari_sp av = avma;
    1384      529236 :   GEN U = mulii(Fp_inv(A,B), A);
    1385      529236 :   return gerepileuptoint(av, Z_chinese_post(a,b,C,U, NULL));
    1386             : }
    1387             : ulong
    1388     2178045 : u_chinese_coprime(ulong a, ulong b, ulong A, ulong B, ulong C)
    1389     2178045 : { return u_chinese_post(a,b,C, A * Fl_inv(A % B,B)); }
    1390             : 
    1391             : /* chinese1 for coprime moduli in Z */
    1392             : static GEN
    1393     2184786 : chinese1_coprime_Z_aux(GEN x, GEN y)
    1394             : {
    1395     2184786 :   GEN z = cgetg(3, t_INTMOD);
    1396     2184786 :   GEN A = gel(x,1), a = gel(x, 2);
    1397     2184786 :   GEN B = gel(y,1), b = gel(y, 2), C = mulii(A,B);
    1398     2184786 :   pari_sp av = avma;
    1399     2184786 :   GEN U = mulii(Fp_inv(A,B), A);
    1400     2184786 :   gel(z,2) = gerepileuptoint(av, Z_chinese_post(a,b,C,U, NULL));
    1401     2184786 :   gel(z,1) = C; return z;
    1402             : }
    1403             : GEN
    1404     2399278 : chinese1_coprime_Z(GEN x) {return gen_chinese(x,chinese1_coprime_Z_aux);}
    1405             : 
    1406             : /*********************************************************************/
    1407             : /**                    MODULAR EXPONENTIATION                       **/
    1408             : /*********************************************************************/
    1409             : /* xa ZV or nv */
    1410             : GEN
    1411     2170985 : ZV_producttree(GEN xa)
    1412             : {
    1413     2170985 :   long n = lg(xa)-1;
    1414     2170985 :   long m = n==1 ? 1: expu(n-1)+1;
    1415     2170986 :   GEN T = cgetg(m+1, t_VEC), t;
    1416             :   long i, j, k;
    1417     2170991 :   t = cgetg(((n+1)>>1)+1, t_VEC);
    1418     2170983 :   if (typ(xa)==t_VECSMALL)
    1419             :   {
    1420     2849274 :     for (j=1, k=1; k<n; j++, k+=2)
    1421     1805005 :       gel(t, j) = muluu(xa[k], xa[k+1]);
    1422     1044269 :     if (k==n) gel(t, j) = utoi(xa[k]);
    1423             :   } else {
    1424     2323895 :     for (j=1, k=1; k<n; j++, k+=2)
    1425     1197198 :       gel(t, j) = mulii(gel(xa,k), gel(xa,k+1));
    1426     1126697 :     if (k==n) gel(t, j) = icopy(gel(xa,k));
    1427             :   }
    1428     2170965 :   gel(T,1) = t;
    1429     3451700 :   for (i=2; i<=m; i++)
    1430             :   {
    1431     1280725 :     GEN u = gel(T, i-1);
    1432     1280725 :     long n = lg(u)-1;
    1433     1280725 :     t = cgetg(((n+1)>>1)+1, t_VEC);
    1434     2857149 :     for (j=1, k=1; k<n; j++, k+=2)
    1435     1576414 :       gel(t, j) = mulii(gel(u, k), gel(u, k+1));
    1436     1280735 :     if (k==n) gel(t, j) = gel(u, k);
    1437     1280735 :     gel(T, i) = t;
    1438             :   }
    1439     2170975 :   return T;
    1440             : }
    1441             : 
    1442             : /* return [A mod P[i], i=1..#P], T = ZV_producttree(P) */
    1443             : GEN
    1444    53532992 : Z_ZV_mod_tree(GEN A, GEN P, GEN T)
    1445             : {
    1446             :   long i,j,k;
    1447    53532992 :   long m = lg(T)-1, n = lg(P)-1;
    1448             :   GEN t;
    1449    53532992 :   GEN Tp = cgetg(m+1, t_VEC);
    1450    53485673 :   gel(Tp, m) = mkvec(modii(A, gmael(T,m,1)));
    1451   112802011 :   for (i=m-1; i>=1; i--)
    1452             :   {
    1453    59398998 :     GEN u = gel(T, i);
    1454    59398998 :     GEN v = gel(Tp, i+1);
    1455    59398998 :     long n = lg(u)-1;
    1456    59398998 :     t = cgetg(n+1, t_VEC);
    1457   147930176 :     for (j=1, k=1; k<n; j++, k+=2)
    1458             :     {
    1459    88655007 :       gel(t, k)   = modii(gel(v, j), gel(u, k));
    1460    88672986 :       gel(t, k+1) = modii(gel(v, j), gel(u, k+1));
    1461             :     }
    1462    59275169 :     if (k==n) gel(t, k) = gel(v, j);
    1463    59275169 :     gel(Tp, i) = t;
    1464             :   }
    1465             :   {
    1466    53403013 :     GEN u = gel(T, i+1);
    1467    53403013 :     GEN v = gel(Tp, i+1);
    1468    53403013 :     long l = lg(u)-1;
    1469    53403013 :     if (typ(P)==t_VECSMALL)
    1470             :     {
    1471    51237348 :       GEN R = cgetg(n+1, t_VECSMALL);
    1472   190788161 :       for (j=1, k=1; j<=l; j++, k+=2)
    1473             :       {
    1474   139178585 :         uel(R,k) = umodiu(gel(v, j), P[k]);
    1475   139523503 :         if (k < n)
    1476   111973078 :           uel(R,k+1) = umodiu(gel(v, j), P[k+1]);
    1477             :       }
    1478    51609576 :       return R;
    1479             :     }
    1480             :     else
    1481             :     {
    1482     2165665 :       GEN R = cgetg(n+1, t_VEC);
    1483     5914667 :       for (j=1, k=1; j<=l; j++, k+=2)
    1484             :       {
    1485     3744081 :         gel(R,k) = modii(gel(v, j), gel(P,k));
    1486     3744077 :         if (k < n)
    1487     2999007 :           gel(R,k+1) = modii(gel(v, j), gel(P,k+1));
    1488             :       }
    1489     2170586 :       return R;
    1490             :     }
    1491             :   }
    1492             : }
    1493             : 
    1494             : /* T = ZV_producttree(P), R = ZV_chinesetree(P,T) */
    1495             : GEN
    1496    35292624 : ZV_chinese_tree(GEN A, GEN P, GEN T, GEN R)
    1497             : {
    1498    35292624 :   long m = lg(T)-1, n = lg(A)-1;
    1499             :   long i,j,k;
    1500    35292624 :   GEN Tp = cgetg(m+1, t_VEC);
    1501    35270293 :   GEN M = gel(T, 1);
    1502    35270293 :   GEN t = cgetg(lg(M), t_VEC);
    1503    35241126 :   if (typ(P)==t_VECSMALL)
    1504             :   {
    1505    82120195 :     for (j=1, k=1; k<n; j++, k+=2)
    1506             :     {
    1507    61163012 :       pari_sp av = avma;
    1508    61163012 :       GEN a = mului(A[k], gel(R,k)), b = mului(A[k+1], gel(R,k+1));
    1509    60982336 :       GEN tj = modii(addii(mului(P[k],b), mului(P[k+1],a)), gel(M,j));
    1510    61071895 :       gel(t, j) = gerepileuptoint(av, tj);
    1511             :     }
    1512    20957183 :     if (k==n) gel(t, j) = modii(mului(A[k], gel(R,k)), gel(M, j));
    1513             :   } else
    1514             :   {
    1515    30386276 :     for (j=1, k=1; k<n; j++, k+=2)
    1516             :     {
    1517    16079264 :       pari_sp av = avma;
    1518    16079264 :       GEN a = mulii(gel(A,k), gel(R,k)), b = mulii(gel(A,k+1), gel(R,k+1));
    1519    16077761 :       GEN tj = modii(addii(mulii(gel(P,k),b), mulii(gel(P,k+1),a)), gel(M,j));
    1520    16103639 :       gel(t, j) = gerepileuptoint(av, tj);
    1521             :     }
    1522    14307012 :     if (k==n) gel(t, j) = modii(mulii(gel(A,k), gel(R,k)), gel(M, j));
    1523             :   }
    1524    35260194 :   gel(Tp, 1) = t;
    1525    67481218 :   for (i=2; i<=m; i++)
    1526             :   {
    1527    32207580 :     GEN u = gel(T, i-1), M = gel(T, i);
    1528    32207580 :     GEN t = cgetg(lg(M), t_VEC);
    1529    32256413 :     GEN v = gel(Tp, i-1);
    1530    32256413 :     long n = lg(v)-1;
    1531    88168843 :     for (j=1, k=1; k<n; j++, k+=2)
    1532             :     {
    1533    55947819 :       pari_sp av = avma;
    1534    55890080 :       gel(t, j) = gerepileuptoint(av, modii(addii(mulii(gel(u, k), gel(v, k+1)),
    1535    55947819 :             mulii(gel(u, k+1), gel(v, k))), gel(M, j)));
    1536             :     }
    1537    32221024 :     if (k==n) gel(t, j) = gel(v, k);
    1538    32221024 :     gel(Tp, i) = t;
    1539             :   }
    1540    35273638 :   return gmael(Tp,m,1);
    1541             : }
    1542             : 
    1543             : static GEN
    1544      937784 : ncV_polint_center_tree(GEN vA, GEN P, GEN T, GEN R, GEN m2)
    1545             : {
    1546      937784 :   long i, l = lg(gel(vA,1)), n = lg(P);
    1547      937784 :   GEN mod = gmael(T, lg(T)-1, 1), V = cgetg(l, t_COL);
    1548    30564388 :   for (i=1; i < l; i++)
    1549             :   {
    1550    29626712 :     pari_sp av = avma;
    1551    29626712 :     GEN c, A = cgetg(n, typ(P));
    1552             :     long j;
    1553   183053771 :     for (j=1; j < n; j++) A[j] = mael(vA,j,i);
    1554    29596419 :     c = Fp_center(ZV_chinese_tree(A, P, T, R), mod, m2);
    1555    29616502 :     gel(V,i) = gerepileuptoint(av, c);
    1556             :   }
    1557      937676 :   return V;
    1558             : }
    1559             : 
    1560             : static GEN
    1561      553476 : nxV_polint_center_tree(GEN vA, GEN P, GEN T, GEN R, GEN m2)
    1562             : {
    1563      553476 :   long i, j, l, n = lg(P);
    1564      553476 :   GEN mod = gmael(T, lg(T)-1, 1), V, w;
    1565      553476 :   w = cgetg(n, t_VECSMALL);
    1566     1847779 :   for(j=1; j<n; j++) w[j] = lg(gel(vA,j));
    1567      553473 :   l = vecsmall_max(w);
    1568      553474 :   V = cgetg(l, t_POL);
    1569      553454 :   V[1] = mael(vA,1,1);
    1570     3225304 :   for (i=2; i < l; i++)
    1571             :   {
    1572     2671837 :     pari_sp av = avma;
    1573     2671837 :     GEN c, A = cgetg(n, typ(P));
    1574     2671569 :     if (typ(P)==t_VECSMALL)
    1575     5973746 :       for (j=1; j < n; j++) A[j] = i < w[j] ? mael(vA,j,i): 0;
    1576             :     else
    1577     3809957 :       for (j=1; j < n; j++) gel(A,j) = i < w[j] ? gmael(vA,j,i): gen_0;
    1578     2671569 :     c = Fp_center(ZV_chinese_tree(A, P, T, R), mod, m2);
    1579     2671814 :     gel(V,i) = gerepileuptoint(av, c);
    1580             :   }
    1581      553467 :   return ZX_renormalize(V, l);
    1582             : }
    1583             : 
    1584             : static GEN
    1585        4614 : nxCV_polint_center_tree(GEN vA, GEN P, GEN T, GEN R, GEN m2)
    1586             : {
    1587        4614 :   long i, j, l = lg(gel(vA,1)), n = lg(P);
    1588        4614 :   GEN A = cgetg(n, t_VEC);
    1589        4614 :   GEN V = cgetg(l, t_COL);
    1590       90903 :   for (i=1; i < l; i++)
    1591             :   {
    1592      335111 :     for (j=1; j < n; j++) gel(A,j) = gmael(vA,j,i);
    1593       86289 :     gel(V,i) = nxV_polint_center_tree(A, P, T, R, m2);
    1594             :   }
    1595        4614 :   return V;
    1596             : }
    1597             : 
    1598             : static GEN
    1599      114988 : polint_chinese(GEN worker, GEN mA, GEN P)
    1600             : {
    1601      114988 :   long cnt, pending, n, i, j, l = lg(gel(mA,1));
    1602             :   struct pari_mt pt;
    1603             :   GEN done, va, M, A;
    1604             :   pari_timer ti;
    1605             : 
    1606      114988 :   if (l == 1) return cgetg(1, t_MAT);
    1607       85936 :   cnt = pending = 0;
    1608       85936 :   n = lg(P);
    1609       85936 :   A = cgetg(n, t_VEC);
    1610       85936 :   va = mkvec(A);
    1611       85936 :   M = cgetg(l, t_MAT);
    1612       85936 :   if (DEBUGLEVEL>4) timer_start(&ti);
    1613       85936 :   if (DEBUGLEVEL>5) err_printf("Start parallel Chinese remainder: ");
    1614       85936 :   mt_queue_start_lim(&pt, worker, l-1);
    1615      676057 :   for (i=1; i<l || pending; i++)
    1616             :   {
    1617             :     long workid;
    1618     2469587 :     for(j=1; j < n; j++) gel(A,j) = gmael(mA,j,i);
    1619      590121 :     mt_queue_submit(&pt, i, i<l? va: NULL);
    1620      590121 :     done = mt_queue_get(&pt, &workid, &pending);
    1621      590121 :     if (done)
    1622             :     {
    1623      555808 :       gel(M,workid) = done;
    1624      555808 :       if (DEBUGLEVEL>5) err_printf("%ld%% ",(++cnt)*100/(l-1));
    1625             :     }
    1626             :   }
    1627       85936 :   if (DEBUGLEVEL>5) err_printf("\n");
    1628       85936 :   if (DEBUGLEVEL>4) timer_printf(&ti, "nmV_chinese_center");
    1629       85936 :   mt_queue_end(&pt);
    1630       85936 :   return M;
    1631             : }
    1632             : 
    1633             : GEN
    1634         840 : nxMV_polint_center_tree_worker(GEN vA, GEN T, GEN R, GEN P, GEN m2)
    1635             : {
    1636         840 :   return nxCV_polint_center_tree(vA, P, T, R, m2);
    1637             : }
    1638             : 
    1639             : static GEN
    1640         430 : nxMV_polint_center_tree_seq(GEN vA, GEN P, GEN T, GEN R, GEN m2)
    1641             : {
    1642         430 :   long i, j, l = lg(gel(vA,1)), n = lg(P);
    1643         430 :   GEN A = cgetg(n, t_VEC);
    1644         430 :   GEN V = cgetg(l, t_MAT);
    1645        4204 :   for (i=1; i < l; i++)
    1646             :   {
    1647       15299 :     for (j=1; j < n; j++) gel(A,j) = gmael(vA,j,i);
    1648        3774 :     gel(V,i) = nxCV_polint_center_tree(A, P, T, R, m2);
    1649             :   }
    1650         430 :   return V;
    1651             : }
    1652             : 
    1653             : static GEN
    1654          90 : nxMV_polint_center_tree(GEN mA, GEN P, GEN T, GEN R, GEN m2)
    1655             : {
    1656          90 :   GEN worker = snm_closure(is_entry("_nxMV_polint_worker"), mkvec4(T, R, P, m2));
    1657          90 :   return polint_chinese(worker, mA, P);
    1658             : }
    1659             : 
    1660             : static GEN
    1661       53115 : nmV_polint_center_tree_seq(GEN vA, GEN P, GEN T, GEN R, GEN m2)
    1662             : {
    1663       53115 :   long i, j, l = lg(gel(vA,1)), n = lg(P);
    1664       53115 :   GEN A = cgetg(n, t_VEC);
    1665       53114 :   GEN V = cgetg(l, t_MAT);
    1666      420230 :   for (i=1; i < l; i++)
    1667             :   {
    1668     2137025 :     for (j=1; j < n; j++) gel(A,j) = gmael(vA,j,i);
    1669      367115 :     gel(V,i) = ncV_polint_center_tree(A, P, T, R, m2);
    1670             :   }
    1671       53115 :   return V;
    1672             : }
    1673             : 
    1674             : GEN
    1675      554898 : nmV_polint_center_tree_worker(GEN vA, GEN T, GEN R, GEN P, GEN m2)
    1676             : {
    1677      554898 :   return ncV_polint_center_tree(vA, P, T, R, m2);
    1678             : }
    1679             : 
    1680             : static GEN
    1681      114898 : nmV_polint_center_tree(GEN mA, GEN P, GEN T, GEN R, GEN m2)
    1682             : {
    1683      114898 :   GEN worker = snm_closure(is_entry("_polint_worker"), mkvec4(T, R, P, m2));
    1684      114898 :   return polint_chinese(worker, mA, P);
    1685             : }
    1686             : 
    1687             : /* return [A mod P[i], i=1..#P] */
    1688             : GEN
    1689           0 : Z_ZV_mod(GEN A, GEN P)
    1690             : {
    1691           0 :   pari_sp av = avma;
    1692           0 :   return gerepilecopy(av, Z_ZV_mod_tree(A, P, ZV_producttree(P)));
    1693             : }
    1694             : /* P a t_VECSMALL */
    1695             : GEN
    1696           0 : Z_nv_mod(GEN A, GEN P)
    1697             : {
    1698           0 :   pari_sp av = avma;
    1699           0 :   return gerepileuptoleaf(av, Z_ZV_mod_tree(A, P, ZV_producttree(P)));
    1700             : }
    1701             : /* B a ZX, T = ZV_producttree(P) */
    1702             : GEN
    1703     1263349 : ZX_nv_mod_tree(GEN B, GEN A, GEN T)
    1704             : {
    1705             :   pari_sp av;
    1706     1263349 :   long i, j, l = lg(B), n = lg(A)-1;
    1707     1263349 :   GEN V = cgetg(n+1, t_VEC);
    1708     6210893 :   for (j=1; j <= n; j++)
    1709             :   {
    1710     4947633 :     gel(V, j) = cgetg(l, t_VECSMALL);
    1711     4947560 :     mael(V, j, 1) = B[1]&VARNBITS;
    1712             :   }
    1713     1263260 :   av = avma;
    1714    12136546 :   for (i=2; i < l; i++)
    1715             :   {
    1716    10873695 :     GEN v = Z_ZV_mod_tree(gel(B, i), A, T);
    1717    75290334 :     for (j=1; j <= n; j++)
    1718    64421866 :       mael(V, j, i) = v[j];
    1719    10868468 :     set_avma(av);
    1720             :   }
    1721     6211100 :   for (j=1; j <= n; j++)
    1722     4948200 :     (void) Flx_renormalize(gel(V, j), l);
    1723     1262900 :   return V;
    1724             : }
    1725             : 
    1726             : static GEN
    1727      169745 : to_ZX(GEN a, long v) { return typ(a)==t_INT? scalarpol(a,v): a; }
    1728             : 
    1729             : GEN
    1730       11939 : ZXX_nv_mod_tree(GEN P, GEN xa, GEN T, long w)
    1731             : {
    1732       11939 :   pari_sp av = avma;
    1733       11939 :   long i, j, l = lg(P), n = lg(xa)-1;
    1734       11939 :   GEN V = cgetg(n+1, t_VEC);
    1735       46965 :   for (j=1; j <= n; j++)
    1736             :   {
    1737       35026 :     gel(V, j) = cgetg(l, t_POL);
    1738       35026 :     mael(V, j, 1) = P[1]&VARNBITS;
    1739             :   }
    1740       99692 :   for (i=2; i < l; i++)
    1741             :   {
    1742       87753 :     GEN v = ZX_nv_mod_tree(to_ZX(gel(P, i), w), xa, T);
    1743      362381 :     for (j=1; j <= n; j++)
    1744      274628 :       gmael(V, j, i) = gel(v,j);
    1745             :   }
    1746       46965 :   for (j=1; j <= n; j++)
    1747       35026 :     (void) FlxX_renormalize(gel(V, j), l);
    1748       11939 :   return gerepilecopy(av, V);
    1749             : }
    1750             : 
    1751             : GEN
    1752        4045 : ZXC_nv_mod_tree(GEN C, GEN xa, GEN T, long w)
    1753             : {
    1754        4045 :   pari_sp av = avma;
    1755        4045 :   long i, j, l = lg(C), n = lg(xa)-1;
    1756        4045 :   GEN V = cgetg(n+1, t_VEC);
    1757       16915 :   for (j = 1; j <= n; j++)
    1758       12870 :     gel(V, j) = cgetg(l, t_COL);
    1759       86037 :   for (i = 1; i < l; i++)
    1760             :   {
    1761       81993 :     GEN v = ZX_nv_mod_tree(to_ZX(gel(C, i), w), xa, T);
    1762      360579 :     for (j = 1; j <= n; j++)
    1763      278587 :       gmael(V, j, i) = gel(v,j);
    1764             :   }
    1765        4044 :   return gerepilecopy(av, V);
    1766             : }
    1767             : 
    1768             : GEN
    1769         430 : ZXM_nv_mod_tree(GEN M, GEN xa, GEN T, long w)
    1770             : {
    1771         430 :   pari_sp av = avma;
    1772         430 :   long i, j, l = lg(M), n = lg(xa)-1;
    1773         430 :   GEN V = cgetg(n+1, t_VEC);
    1774        2083 :   for (j=1; j <= n; j++)
    1775        1653 :     gel(V, j) = cgetg(l, t_MAT);
    1776        4204 :   for (i=1; i < l; i++)
    1777             :   {
    1778        3774 :     GEN v = ZXC_nv_mod_tree(gel(M, i), xa, T, w);
    1779       15299 :     for (j=1; j <= n; j++)
    1780       11525 :       gmael(V, j, i) = gel(v,j);
    1781             :   }
    1782         430 :   return gerepilecopy(av, V);
    1783             : }
    1784             : 
    1785             : GEN
    1786      941093 : ZV_nv_mod_tree(GEN B, GEN A, GEN T)
    1787             : {
    1788             :   pari_sp av;
    1789      941093 :   long i, j, l = lg(B), n = lg(A)-1;
    1790      941093 :   GEN V = cgetg(n+1, t_VEC);
    1791     4667977 :   for (j=1; j <= n; j++) gel(V, j) = cgetg(l, t_VECSMALL);
    1792      941028 :   av = avma;
    1793    41408892 :   for (i=1; i < l; i++)
    1794             :   {
    1795    40476190 :     GEN v = Z_ZV_mod_tree(gel(B, i), A, T);
    1796   228740005 :     for (j=1; j <= n; j++) mael(V, j, i) = v[j];
    1797    40440222 :     set_avma(av);
    1798             :   }
    1799      932702 :   return V;
    1800             : }
    1801             : 
    1802             : GEN
    1803       80998 : ZM_nv_mod_tree(GEN M, GEN xa, GEN T)
    1804             : {
    1805       80998 :   pari_sp av = avma;
    1806       80998 :   long i, j, l = lg(M), n = lg(xa)-1;
    1807       80998 :   GEN V = cgetg(n+1, t_VEC);
    1808      317049 :   for (j=1; j <= n; j++) gel(V, j) = cgetg(l, t_MAT);
    1809     1021856 :   for (i=1; i < l; i++)
    1810             :   {
    1811      940866 :     GEN v = ZV_nv_mod_tree(gel(M, i), xa, T);
    1812     4667441 :     for (j=1; j <= n; j++) gmael(V, j, i) = gel(v,j);
    1813             :   }
    1814       80990 :   return gerepilecopy(av, V);
    1815             : }
    1816             : 
    1817             : static GEN
    1818     2166733 : ZV_sqr(GEN z)
    1819             : {
    1820     2166733 :   long i,l = lg(z);
    1821     2166733 :   GEN x = cgetg(l, t_VEC);
    1822     2166753 :   if (typ(z)==t_VECSMALL)
    1823     5041388 :     for (i=1; i<l; i++) gel(x,i) = sqru(z[i]);
    1824             :   else
    1825     3842661 :     for (i=1; i<l; i++) gel(x,i) = sqri(gel(z,i));
    1826     2166709 :   return x;
    1827             : }
    1828             : 
    1829             : static GEN
    1830    11101186 : ZT_sqr(GEN x)
    1831             : {
    1832    11101186 :   if (typ(x) == t_INT) return sqri(x);
    1833    14542189 :   pari_APPLY_type(t_VEC, ZT_sqr(gel(x,i)))
    1834             : }
    1835             : 
    1836             : static GEN
    1837     2166730 : ZV_invdivexact(GEN y, GEN x)
    1838             : {
    1839     2166730 :   long i, l = lg(y);
    1840     2166730 :   GEN z = cgetg(l,t_VEC);
    1841     2166766 :   if (typ(x)==t_VECSMALL)
    1842     5041117 :     for (i=1; i<l; i++)
    1843             :     {
    1844     3997236 :       pari_sp av = avma;
    1845     3997236 :       ulong a = Fl_inv(umodiu(diviuexact(gel(y,i),x[i]), x[i]), x[i]);
    1846     3997493 :       set_avma(av); gel(z,i) = utoi(a);
    1847             :     }
    1848             :   else
    1849     3842774 :     for (i=1; i<l; i++)
    1850     2719903 :       gel(z,i) = Fp_inv(diviiexact(gel(y,i), gel(x,i)), gel(x,i));
    1851     2166752 :   return z;
    1852             : }
    1853             : 
    1854             : /* P t_VECSMALL or t_VEC of t_INT  */
    1855             : GEN
    1856     2166741 : ZV_chinesetree(GEN P, GEN T)
    1857             : {
    1858     2166741 :   GEN T2 = ZT_sqr(T), P2 = ZV_sqr(P);
    1859     2166729 :   GEN mod = gmael(T,lg(T)-1,1);
    1860     2166729 :   return ZV_invdivexact(Z_ZV_mod_tree(mod, P2, T2), P);
    1861             : }
    1862             : 
    1863             : static GEN
    1864      646491 : gc_chinese(pari_sp av, GEN T, GEN a, GEN *pt_mod)
    1865             : {
    1866      646491 :   if (!pt_mod)
    1867       12209 :     return gerepileupto(av, a);
    1868             :   else
    1869             :   {
    1870      634282 :     GEN mod = gmael(T, lg(T)-1, 1);
    1871      634282 :     gerepileall(av, 2, &a, &mod);
    1872      634281 :     *pt_mod = mod;
    1873      634281 :     return a;
    1874             :   }
    1875             : }
    1876             : 
    1877             : GEN
    1878      162125 : ZV_chinese_center(GEN A, GEN P, GEN *pt_mod)
    1879             : {
    1880      162125 :   pari_sp av = avma;
    1881      162125 :   GEN T = ZV_producttree(P);
    1882      162126 :   GEN R = ZV_chinesetree(P, T);
    1883      162126 :   GEN a = ZV_chinese_tree(A, P, T, R);
    1884      162126 :   GEN mod = gmael(T, lg(T)-1, 1);
    1885      162126 :   GEN ca = Fp_center(a, mod, shifti(mod,-1));
    1886      162126 :   return gc_chinese(av, T, ca, pt_mod);
    1887             : }
    1888             : 
    1889             : GEN
    1890        5095 : ZV_chinese(GEN A, GEN P, GEN *pt_mod)
    1891             : {
    1892        5095 :   pari_sp av = avma;
    1893        5095 :   GEN T = ZV_producttree(P);
    1894        5095 :   GEN R = ZV_chinesetree(P, T);
    1895        5095 :   GEN a = ZV_chinese_tree(A, P, T, R);
    1896        5095 :   return gc_chinese(av, T, a, pt_mod);
    1897             : }
    1898             : 
    1899             : GEN
    1900      113129 : nxV_chinese_center_tree(GEN A, GEN P, GEN T, GEN R)
    1901             : {
    1902      113129 :   pari_sp av = avma;
    1903      113129 :   GEN m2 = shifti(gmael(T, lg(T)-1, 1), -1);
    1904      113129 :   GEN a = nxV_polint_center_tree(A, P, T, R, m2);
    1905      113129 :   return gerepileupto(av, a);
    1906             : }
    1907             : 
    1908             : GEN
    1909      354046 : nxV_chinese_center(GEN A, GEN P, GEN *pt_mod)
    1910             : {
    1911      354046 :   pari_sp av = avma;
    1912      354046 :   GEN T = ZV_producttree(P);
    1913      354046 :   GEN R = ZV_chinesetree(P, T);
    1914      354045 :   GEN m2 = shifti(gmael(T, lg(T)-1, 1), -1);
    1915      354044 :   GEN a = nxV_polint_center_tree(A, P, T, R, m2);
    1916      354045 :   return gc_chinese(av, T, a, pt_mod);
    1917             : }
    1918             : 
    1919             : GEN
    1920       10237 : ncV_chinese_center(GEN A, GEN P, GEN *pt_mod)
    1921             : {
    1922       10237 :   pari_sp av = avma;
    1923       10237 :   GEN T = ZV_producttree(P);
    1924       10237 :   GEN R = ZV_chinesetree(P, T);
    1925       10237 :   GEN m2 = shifti(gmael(T, lg(T)-1, 1), -1);
    1926       10237 :   GEN a = ncV_polint_center_tree(A, P, T, R, m2);
    1927       10237 :   return gc_chinese(av, T, a, pt_mod);
    1928             : }
    1929             : 
    1930             : GEN
    1931        5530 : ncV_chinese_center_tree(GEN A, GEN P, GEN T, GEN R)
    1932             : {
    1933        5530 :   pari_sp av = avma;
    1934        5530 :   GEN m2 = shifti(gmael(T, lg(T)-1, 1), -1);
    1935        5530 :   GEN a = ncV_polint_center_tree(A, P, T, R, m2);
    1936        5530 :   return gerepileupto(av, a);
    1937             : }
    1938             : 
    1939             : GEN
    1940           0 : nmV_chinese_center_tree(GEN A, GEN P, GEN T, GEN R)
    1941             : {
    1942           0 :   pari_sp av = avma;
    1943           0 :   GEN m2 = shifti(gmael(T, lg(T)-1, 1), -1);
    1944           0 :   GEN a = nmV_polint_center_tree(A, P, T, R, m2);
    1945           0 :   return gerepileupto(av, a);
    1946             : }
    1947             : 
    1948             : GEN
    1949       53115 : nmV_chinese_center_tree_seq(GEN A, GEN P, GEN T, GEN R)
    1950             : {
    1951       53115 :   pari_sp av = avma;
    1952       53115 :   GEN m2 = shifti(gmael(T, lg(T)-1, 1), -1);
    1953       53115 :   GEN a = nmV_polint_center_tree_seq(A, P, T, R, m2);
    1954       53115 :   return gerepileupto(av, a);
    1955             : }
    1956             : 
    1957             : GEN
    1958      114898 : nmV_chinese_center(GEN A, GEN P, GEN *pt_mod)
    1959             : {
    1960      114898 :   pari_sp av = avma;
    1961      114898 :   GEN T = ZV_producttree(P);
    1962      114898 :   GEN R = ZV_chinesetree(P, T);
    1963      114898 :   GEN m2 = shifti(gmael(T, lg(T)-1, 1), -1);
    1964      114898 :   GEN a = nmV_polint_center_tree(A, P, T, R, m2);
    1965      114898 :   return gc_chinese(av, T, a, pt_mod);
    1966             : }
    1967             : 
    1968             : GEN
    1969           0 : nxCV_chinese_center_tree(GEN A, GEN P, GEN T, GEN R)
    1970             : {
    1971           0 :   pari_sp av = avma;
    1972           0 :   GEN m2 = shifti(gmael(T, lg(T)-1, 1), -1);
    1973           0 :   GEN a = nxCV_polint_center_tree(A, P, T, R, m2);
    1974           0 :   return gerepileupto(av, a);
    1975             : }
    1976             : 
    1977             : GEN
    1978           0 : nxCV_chinese_center(GEN A, GEN P, GEN *pt_mod)
    1979             : {
    1980           0 :   pari_sp av = avma;
    1981           0 :   GEN T = ZV_producttree(P);
    1982           0 :   GEN R = ZV_chinesetree(P, T);
    1983           0 :   GEN m2 = shifti(gmael(T, lg(T)-1, 1), -1);
    1984           0 :   GEN a = nxCV_polint_center_tree(A, P, T, R, m2);
    1985           0 :   return gc_chinese(av, T, a, pt_mod);
    1986             : }
    1987             : 
    1988             : GEN
    1989         430 : nxMV_chinese_center_tree_seq(GEN A, GEN P, GEN T, GEN R)
    1990             : {
    1991         430 :   pari_sp av = avma;
    1992         430 :   GEN m2 = shifti(gmael(T, lg(T)-1, 1), -1);
    1993         430 :   GEN a = nxMV_polint_center_tree_seq(A, P, T, R, m2);
    1994         430 :   return gerepileupto(av, a);
    1995             : }
    1996             : 
    1997             : GEN
    1998          90 : nxMV_chinese_center(GEN A, GEN P, GEN *pt_mod)
    1999             : {
    2000          90 :   pari_sp av = avma;
    2001          90 :   GEN T = ZV_producttree(P);
    2002          90 :   GEN R = ZV_chinesetree(P, T);
    2003          90 :   GEN m2 = shifti(gmael(T, lg(T)-1, 1), -1);
    2004          90 :   GEN a = nxMV_polint_center_tree(A, P, T, R, m2);
    2005          90 :   return gc_chinese(av, T, a, pt_mod);
    2006             : }
    2007             : 
    2008             : /**********************************************************************
    2009             :  **                    Powering  over (Z/NZ)^*, small N              **
    2010             :  **********************************************************************/
    2011             : 
    2012             : /* 2^n mod p; assume n > 1 */
    2013             : static ulong
    2014    10406985 : Fl_2powu_pre(ulong n, ulong p, ulong pi)
    2015             : {
    2016    10406985 :   ulong y = 2;
    2017    10406985 :   int j = 1+bfffo(n);
    2018             :   /* normalize, i.e set highest bit to 1 (we know n != 0) */
    2019    10406985 :   n<<=j; j = BITS_IN_LONG-j; /* first bit is now implicit */
    2020   444967317 :   for (; j; n<<=1,j--)
    2021             :   {
    2022   434611911 :     y = Fl_sqr_pre(y,p,pi);
    2023   434675646 :     if (n & HIGHBIT) y = Fl_double(y, p);
    2024             :   }
    2025    10355406 :   return y;
    2026             : }
    2027             : 
    2028             : /* 2^n mod p; assume n > 1 and !(p & HIGHMASK) */
    2029             : static ulong
    2030     4172176 : Fl_2powu(ulong n, ulong p)
    2031             : {
    2032     4172176 :   ulong y = 2;
    2033     4172176 :   int j = 1+bfffo(n);
    2034             :   /* normalize, i.e set highest bit to 1 (we know n != 0) */
    2035     4172176 :   n<<=j; j = BITS_IN_LONG-j; /* first bit is now implicit */
    2036    23878999 :   for (; j; n<<=1,j--)
    2037             :   {
    2038    19706820 :     y = (y*y) % p;
    2039    19706820 :     if (n & HIGHBIT) y = Fl_double(y, p);
    2040             :   }
    2041     4172179 :   return y;
    2042             : }
    2043             : 
    2044             : /* allow pi = 0 */
    2045             : ulong
    2046    92138713 : Fl_powu_pre(ulong x, ulong n0, ulong p, ulong pi)
    2047             : {
    2048             :   ulong y, z, n;
    2049    92138713 :   if (!pi) return Fl_powu(x, n0, p);
    2050    89806538 :   if (n0 <= 1)
    2051             :   { /* frequent special cases */
    2052     5520693 :     if (n0 == 1) return x;
    2053       81176 :     if (n0 == 0) return 1;
    2054             :   }
    2055    84285847 :   if (x <= 2)
    2056             :   {
    2057    10450806 :     if (x == 2) return Fl_2powu_pre(n0, p, pi);
    2058       42689 :     return x; /* 0 or 1 */
    2059             :   }
    2060    73835041 :   y = 1; z = x; n = n0;
    2061             :   for(;;)
    2062             :   {
    2063   529918699 :     if (n&1) y = Fl_mul_pre(y,z,p,pi);
    2064   529808295 :     n>>=1; if (!n) return y;
    2065   455835890 :     z = Fl_sqr_pre(z,p,pi);
    2066             :   }
    2067             : }
    2068             : 
    2069             : ulong
    2070   144946162 : Fl_powu(ulong x, ulong n0, ulong p)
    2071             : {
    2072             :   ulong y, z, n;
    2073   144946162 :   if (n0 <= 2)
    2074             :   { /* frequent special cases */
    2075    65326515 :     if (n0 == 2) return Fl_sqr(x,p);
    2076    31910467 :     if (n0 == 1) return x;
    2077     1883736 :     if (n0 == 0) return 1;
    2078             :   }
    2079    79588034 :   if (x <= 1) return x; /* 0 or 1 */
    2080    79167748 :   if (p & HIGHMASK)
    2081     7881972 :     return Fl_powu_pre(x, n0, p, get_Fl_red(p));
    2082    71285776 :   if (x == 2) return Fl_2powu(n0, p);
    2083    67113591 :   y = 1; z = x; n = n0;
    2084             :   for(;;)
    2085             :   {
    2086   354463454 :     if (n&1) y = (y*z) % p;
    2087   354463454 :     n>>=1; if (!n) return y;
    2088   287349863 :     z = (z*z) % p;
    2089             :   }
    2090             : }
    2091             : 
    2092             : /* Reduce data dependency to maximize internal parallelism; allow pi = 0 */
    2093             : GEN
    2094    12768563 : Fl_powers_pre(ulong x, long n, ulong p, ulong pi)
    2095             : {
    2096             :   long i, k;
    2097    12768563 :   GEN z = cgetg(n + 2, t_VECSMALL);
    2098    12757698 :   z[1] = 1; if (n == 0) return z;
    2099    12757698 :   z[2] = x;
    2100    12757698 :   if (pi)
    2101             :   {
    2102    90048449 :     for (i = 3, k=2; i <= n; i+=2, k++)
    2103             :     {
    2104    77494034 :       z[i] = Fl_sqr_pre(z[k], p, pi);
    2105    77499431 :       z[i+1] = Fl_mul_pre(z[k], z[k+1], p, pi);
    2106             :     }
    2107    12554415 :     if (i==n+1) z[i] = Fl_sqr_pre(z[k], p, pi);
    2108             :   }
    2109      212378 :   else if (p & HIGHMASK)
    2110             :   {
    2111           0 :     for (i = 3, k=2; i <= n; i+=2, k++)
    2112             :     {
    2113           0 :       z[i] = Fl_sqr(z[k], p);
    2114           0 :       z[i+1] = Fl_mul(z[k], z[k+1], p);
    2115             :     }
    2116           0 :     if (i==n+1) z[i] = Fl_sqr(z[k], p);
    2117             :   }
    2118             :   else
    2119   400517421 :     for (i = 2; i <= n; i++) z[i+1] = (z[i] * x) % p;
    2120    12768124 :   return z;
    2121             : }
    2122             : 
    2123             : GEN
    2124      295824 : Fl_powers(ulong x, long n, ulong p)
    2125             : {
    2126      295824 :   return Fl_powers_pre(x, n, p, (p & HIGHMASK)? get_Fl_red(p): 0);
    2127             : }
    2128             : 
    2129             : /**********************************************************************
    2130             :  **                    Powering  over (Z/NZ)^*, large N              **
    2131             :  **********************************************************************/
    2132             : 
    2133             : static GEN
    2134     4565239 : Fp_dblsqr(GEN x, GEN N)
    2135             : {
    2136     4565239 :   GEN z = shifti(Fp_sqr(x, N), 1);
    2137     4466924 :   return cmpii(z, N) >= 0? subii(z, N): z;
    2138             : }
    2139             : 
    2140             : typedef struct muldata {
    2141             :   GEN (*sqr)(void * E, GEN x);
    2142             :   GEN (*mul)(void * E, GEN x, GEN y);
    2143             :   GEN (*mul2)(void * E, GEN x);
    2144             : } muldata;
    2145             : 
    2146             : /* modified Barrett reduction with one fold */
    2147             : /* See Fast Modular Reduction, W. Hasenplaugh, G. Gaubatz, V. Gopal, ARITH 18 */
    2148             : 
    2149             : static GEN
    2150       12951 : Fp_invmBarrett(GEN p, long s)
    2151             : {
    2152       12951 :   GEN R, Q = dvmdii(int2n(3*s),p,&R);
    2153       12951 :   return mkvec2(Q,R);
    2154             : }
    2155             : 
    2156             : /* a <= (N-1)^2, 2^(2s-2) <= N < 2^(2s). Return 0 <= r < N such that
    2157             :  * a = r (mod N) */
    2158             : static GEN
    2159     6401240 : Fp_rem_mBarrett(GEN a, GEN B, long s, GEN N)
    2160             : {
    2161     6401240 :   pari_sp av = avma;
    2162     6401240 :   GEN P = gel(B, 1), Q = gel(B, 2); /* 2^(3s) = P N + Q, 0 <= Q < N */
    2163     6401240 :   long t = expi(P)+1; /* 2^(t-1) <= P < 2^t */
    2164     6401240 :   GEN u = shifti(a, -3*s), v = remi2n(a, 3*s); /* a = 2^(3s)u + v */
    2165     6401240 :   GEN A = addii(v, mulii(Q,u)); /* 0 <= A < 2^(3s+1) */
    2166     6401240 :   GEN q = shifti(mulii(shifti(A, t-3*s), P), -t); /* A/N - 4 < q <= A/N */
    2167     6401240 :   GEN r = subii(A, mulii(q, N));
    2168     6401240 :   GEN sr= subii(r,N);     /* 0 <= r < 4*N */
    2169     6401240 :   if (signe(sr)<0) return gerepileuptoint(av, r);
    2170     3343618 :   r=sr; sr = subii(r,N);  /* 0 <= r < 3*N */
    2171     3343618 :   if (signe(sr)<0) return gerepileuptoint(av, r);
    2172      105378 :   r=sr; sr = subii(r,N);  /* 0 <= r < 2*N */
    2173      105378 :   return gerepileuptoint(av, signe(sr)>=0 ? sr:r);
    2174             : }
    2175             : 
    2176             : /* Montgomery reduction */
    2177             : 
    2178             : INLINE ulong
    2179      864342 : init_montdata(GEN N) { return (ulong) -invmod2BIL(mod2BIL(N)); }
    2180             : 
    2181             : struct montred
    2182             : {
    2183             :   GEN N;
    2184             :   ulong inv;
    2185             : };
    2186             : 
    2187             : /* Montgomery reduction */
    2188             : static GEN
    2189    61131223 : _sqr_montred(void * E, GEN x)
    2190             : {
    2191    61131223 :   struct montred * D = (struct montred *) E;
    2192    61131223 :   return red_montgomery(sqri(x), D->N, D->inv);
    2193             : }
    2194             : 
    2195             : /* Montgomery reduction */
    2196             : static GEN
    2197     7371143 : _mul_montred(void * E, GEN x, GEN y)
    2198             : {
    2199     7371143 :   struct montred * D = (struct montred *) E;
    2200     7371143 :   return red_montgomery(mulii(x, y), D->N, D->inv);
    2201             : }
    2202             : 
    2203             : static GEN
    2204     6863019 : _mul2_montred(void * E, GEN x)
    2205             : {
    2206     6863019 :   struct montred * D = (struct montred *) E;
    2207     6863019 :   GEN z = shifti(_sqr_montred(E, x), 1);
    2208     6862567 :   long l = lgefint(D->N);
    2209     7289776 :   while (lgefint(z) > l) z = subii(z, D->N);
    2210     6862808 :   return z;
    2211             : }
    2212             : 
    2213             : static GEN
    2214    20454433 : _sqr_remii(void* N, GEN x)
    2215    20454433 : { return remii(sqri(x), (GEN) N); }
    2216             : 
    2217             : static GEN
    2218     1598651 : _mul_remii(void* N, GEN x, GEN y)
    2219     1598651 : { return remii(mulii(x, y), (GEN) N); }
    2220             : 
    2221             : static GEN
    2222     3221120 : _mul2_remii(void* N, GEN x)
    2223     3221120 : { return Fp_dblsqr(x, (GEN) N); }
    2224             : 
    2225             : struct redbarrett
    2226             : {
    2227             :   GEN iM, N;
    2228             :   long s;
    2229             : };
    2230             : 
    2231             : static GEN
    2232     5626446 : _sqr_remiibar(void *E, GEN x)
    2233             : {
    2234     5626446 :   struct redbarrett * D = (struct redbarrett *) E;
    2235     5626446 :   return Fp_rem_mBarrett(sqri(x), D->iM, D->s, D->N);
    2236             : }
    2237             : 
    2238             : static GEN
    2239      774794 : _mul_remiibar(void *E, GEN x, GEN y)
    2240             : {
    2241      774794 :   struct redbarrett * D = (struct redbarrett *) E;
    2242      774794 :   return Fp_rem_mBarrett(mulii(x, y), D->iM, D->s, D->N);
    2243             : }
    2244             : 
    2245             : static GEN
    2246     1342530 : _mul2_remiibar(void *E, GEN x)
    2247             : {
    2248     1342530 :   struct redbarrett * D = (struct redbarrett *) E;
    2249     1342530 :   return Fp_dblsqr(x, D->N);
    2250             : }
    2251             : 
    2252             : static long
    2253     1084957 : Fp_select_red(GEN *y, ulong k, GEN N, long lN, muldata *D, void **pt_E)
    2254             : {
    2255     1084957 :   if (lN >= Fp_POW_BARRETT_LIMIT && (k==0 || ((double)k)*expi(*y) > 2 + expi(N)))
    2256             :   {
    2257       12951 :     struct redbarrett * E = (struct redbarrett *) stack_malloc(sizeof(struct redbarrett));
    2258       12951 :     D->sqr = &_sqr_remiibar;
    2259       12951 :     D->mul = &_mul_remiibar;
    2260       12951 :     D->mul2 = &_mul2_remiibar;
    2261       12951 :     E->N = N;
    2262       12951 :     E->s = 1+(expi(N)>>1);
    2263       12951 :     E->iM = Fp_invmBarrett(N, E->s);
    2264       12951 :     *pt_E = (void*) E;
    2265       12951 :     return 0;
    2266             :   }
    2267     1072006 :   else if (mod2(N) && lN < Fp_POW_REDC_LIMIT)
    2268             :   {
    2269      864341 :     struct montred * E = (struct montred *) stack_malloc(sizeof(struct montred));
    2270      864341 :     *y = remii(shifti(*y, bit_accuracy(lN)), N);
    2271      864345 :     D->sqr = &_sqr_montred;
    2272      864345 :     D->mul = &_mul_montred;
    2273      864345 :     D->mul2 = &_mul2_montred;
    2274      864345 :     E->N = N;
    2275      864345 :     E->inv = init_montdata(N);
    2276      864335 :     *pt_E = (void*) E;
    2277      864335 :     return 1;
    2278             :   }
    2279             :   else
    2280             :   {
    2281      207660 :     D->sqr = &_sqr_remii;
    2282      207660 :     D->mul = &_mul_remii;
    2283      207660 :     D->mul2 = &_mul2_remii;
    2284      207660 :     *pt_E = (void*) N;
    2285      207660 :     return 0;
    2286             :   }
    2287             : }
    2288             : 
    2289             : GEN
    2290     3107374 : Fp_powu(GEN A, ulong k, GEN N)
    2291             : {
    2292     3107374 :   long lN = lgefint(N);
    2293             :   int base_is_2, use_montgomery;
    2294             :   muldata D;
    2295             :   void *E;
    2296             :   pari_sp av;
    2297             : 
    2298     3107374 :   if (lN == 3) {
    2299     1443355 :     ulong n = uel(N,2);
    2300     1443355 :     return utoi( Fl_powu(umodiu(A, n), k, n) );
    2301             :   }
    2302     1664019 :   if (k <= 2)
    2303             :   { /* frequent special cases */
    2304      837425 :     if (k == 2) return Fp_sqr(A,N);
    2305      279860 :     if (k == 1) return A;
    2306           0 :     if (k == 0) return gen_1;
    2307             :   }
    2308      826593 :   av = avma; A = modii(A,N);
    2309      826594 :   base_is_2 = 0;
    2310      826594 :   if (lgefint(A) == 3) switch(A[2])
    2311             :   {
    2312         832 :     case 1: set_avma(av); return gen_1;
    2313       42003 :     case 2:  base_is_2 = 1; break;
    2314             :   }
    2315             : 
    2316             :   /* TODO: Move this out of here and use for general modular computations */
    2317      825762 :   use_montgomery = Fp_select_red(&A, k, N, lN, &D, &E);
    2318      825762 :   if (base_is_2)
    2319       42003 :     A = gen_powu_fold_i(A, k, E, D.sqr, D.mul2);
    2320             :   else
    2321      783759 :     A = gen_powu_i(A, k, E, D.sqr, D.mul);
    2322      825762 :   if (use_montgomery)
    2323             :   {
    2324      704844 :     A = red_montgomery(A, N, ((struct montred *) E)->inv);
    2325      704844 :     if (cmpii(A, N) >= 0) A = subii(A,N);
    2326             :   }
    2327      825762 :   return gerepileuptoint(av, A);
    2328             : }
    2329             : 
    2330             : GEN
    2331       29410 : Fp_pows(GEN A, long k, GEN N)
    2332             : {
    2333       29410 :   if (lgefint(N) == 3) {
    2334       14921 :     ulong n = N[2];
    2335       14921 :     ulong a = umodiu(A, n);
    2336       14921 :     if (k < 0) {
    2337         133 :       a = Fl_inv(a, n);
    2338         133 :       k = -k;
    2339             :     }
    2340       14921 :     return utoi( Fl_powu(a, (ulong)k, n) );
    2341             :   }
    2342       14489 :   if (k < 0) { A = Fp_inv(A, N); k = -k; };
    2343       14489 :   return Fp_powu(A, (ulong)k, N);
    2344             : }
    2345             : 
    2346             : /* A^K mod N */
    2347             : GEN
    2348    35964868 : Fp_pow(GEN A, GEN K, GEN N)
    2349             : {
    2350             :   pari_sp av;
    2351    35964868 :   long s, lN = lgefint(N), sA, sy;
    2352             :   int base_is_2, use_montgomery;
    2353             :   GEN y;
    2354             :   muldata D;
    2355             :   void *E;
    2356             : 
    2357    35964868 :   s = signe(K);
    2358    35964868 :   if (!s) return dvdii(A,N)? gen_0: gen_1;
    2359    34933592 :   if (lN == 3 && lgefint(K) == 3)
    2360             :   {
    2361    33979222 :     ulong n = N[2], a = umodiu(A, n);
    2362    33979484 :     if (s < 0) a = Fl_inv(a, n);
    2363    33979548 :     if (a <= 1) return utoi(a); /* 0 or 1 */
    2364    30481069 :     return utoi(Fl_powu(a, uel(K,2), n));
    2365             :   }
    2366             : 
    2367      954370 :   av = avma;
    2368      954370 :   if (s < 0) y = Fp_inv(A,N);
    2369             :   else
    2370             :   {
    2371      952415 :     y = modii(A,N);
    2372      952616 :     if (!signe(y)) { set_avma(av); return gen_0; }
    2373             :   }
    2374      954571 :   if (lgefint(K) == 3) return gerepileuptoint(av, Fp_powu(y, K[2], N));
    2375             : 
    2376      259365 :   base_is_2 = 0;
    2377      259365 :   sy = abscmpii(y, shifti(N,-1)) > 0;
    2378      259387 :   if (sy) y = subii(N,y);
    2379      259398 :   sA = sy && mod2(K);
    2380      259398 :   if (lgefint(y) == 3) switch(y[2])
    2381             :   {
    2382         200 :     case 1:  set_avma(av); return sA ? subis(N,1): gen_1;
    2383      132923 :     case 2:  base_is_2 = 1; break;
    2384             :   }
    2385             : 
    2386             :   /* TODO: Move this out of here and use for general modular computations */
    2387      259198 :   use_montgomery = Fp_select_red(&y, 0UL, N, lN, &D, &E);
    2388      259185 :   if (base_is_2)
    2389      132923 :     y = gen_pow_fold_i(y, K, E, D.sqr, D.mul2);
    2390             :   else
    2391      126262 :     y = gen_pow_i(y, K, E, D.sqr, D.mul);
    2392      259205 :   if (use_montgomery)
    2393             :   {
    2394      159509 :     y = red_montgomery(y, N, ((struct montred *) E)->inv);
    2395      159508 :     if (cmpii(y,N) >= 0) y = subii(y,N);
    2396             :   }
    2397      259203 :   if (sA) y = subii(N, y);
    2398      259203 :   return gerepileuptoint(av,y);
    2399             : }
    2400             : 
    2401             : static GEN
    2402    14207608 : _Fp_mul(void *E, GEN x, GEN y) { return Fp_mul(x,y,(GEN)E); }
    2403             : static GEN
    2404     8134365 : _Fp_sqr(void *E, GEN x) { return Fp_sqr(x,(GEN)E); }
    2405             : static GEN
    2406       57564 : _Fp_one(void *E) { (void) E; return gen_1; }
    2407             : 
    2408             : GEN
    2409          77 : Fp_pow_init(GEN x, GEN n, long k, GEN p)
    2410          77 : { return gen_pow_init(x, n, k, (void*)p, &_Fp_sqr, &_Fp_mul); }
    2411             : 
    2412             : GEN
    2413       54096 : Fp_pow_table(GEN R, GEN n, GEN p)
    2414       54096 : { return gen_pow_table(R, n, (void*)p, &_Fp_one, &_Fp_mul); }
    2415             : 
    2416             : GEN
    2417        5924 : Fp_powers(GEN x, long n, GEN p)
    2418             : {
    2419        5924 :   if (lgefint(p) == 3)
    2420        2456 :     return Flv_to_ZV(Fl_powers(umodiu(x, uel(p, 2)), n, uel(p, 2)));
    2421        3468 :   return gen_powers(x, n, 1, (void*)p, _Fp_sqr, _Fp_mul, _Fp_one);
    2422             : }
    2423             : 
    2424             : GEN
    2425         497 : FpV_prod(GEN V, GEN p) { return gen_product(V, (void *)p, &_Fp_mul); }
    2426             : 
    2427             : static GEN
    2428    27743391 : _Fp_pow(void *E, GEN x, GEN n) { return Fp_pow(x,n,(GEN)E); }
    2429             : static GEN
    2430         154 : _Fp_rand(void *E) { return addiu(randomi(subiu((GEN)E,1)),1); }
    2431             : 
    2432             : static GEN Fp_easylog(void *E, GEN a, GEN g, GEN ord);
    2433             : static const struct bb_group Fp_star={_Fp_mul,_Fp_pow,_Fp_rand,hash_GEN,
    2434             :                                       equalii,equali1,Fp_easylog};
    2435             : 
    2436             : static GEN
    2437      782379 : _Fp_red(void *E, GEN x) { return Fp_red(x, (GEN)E); }
    2438             : static GEN
    2439      936403 : _Fp_add(void *E, GEN x, GEN y) { (void) E; return addii(x,y); }
    2440             : static GEN
    2441      844296 : _Fp_neg(void *E, GEN x) { (void) E; return negi(x); }
    2442             : static GEN
    2443      513167 : _Fp_rmul(void *E, GEN x, GEN y) { (void) E; return mulii(x,y); }
    2444             : static GEN
    2445       33113 : _Fp_inv(void *E, GEN x) { return Fp_inv(x,(GEN)E); }
    2446             : static int
    2447      222261 : _Fp_equal0(GEN x) { return signe(x)==0; }
    2448             : static GEN
    2449       28726 : _Fp_s(void *E, long x) { (void) E; return stoi(x); }
    2450             : 
    2451             : static const struct bb_field Fp_field={_Fp_red,_Fp_add,_Fp_rmul,_Fp_neg,
    2452             :                                         _Fp_inv,_Fp_equal0,_Fp_s};
    2453             : 
    2454        7093 : const struct bb_field *get_Fp_field(void **E, GEN p)
    2455        7093 : { *E = (void*)p; return &Fp_field; }
    2456             : 
    2457             : /*********************************************************************/
    2458             : /**               ORDER of INTEGERMOD x  in  (Z/nZ)*                **/
    2459             : /*********************************************************************/
    2460             : ulong
    2461      481619 : Fl_order(ulong a, ulong o, ulong p)
    2462             : {
    2463      481619 :   pari_sp av = avma;
    2464             :   GEN m, P, E;
    2465             :   long i;
    2466      481619 :   if (a==1) return 1;
    2467      438733 :   if (!o) o = p-1;
    2468      438733 :   m = factoru(o);
    2469      438733 :   P = gel(m,1);
    2470      438733 :   E = gel(m,2);
    2471     1250931 :   for (i = lg(P)-1; i; i--)
    2472             :   {
    2473      812198 :     ulong j, l = P[i], e = E[i], t = o / upowuu(l,e), y = Fl_powu(a, t, p);
    2474      812198 :     if (y == 1) o = t;
    2475      772762 :     else for (j = 1; j < e; j++)
    2476             :     {
    2477      385554 :       y = Fl_powu(y, l, p);
    2478      385554 :       if (y == 1) { o = t *  upowuu(l, j); break; }
    2479             :     }
    2480             :   }
    2481      438733 :   return gc_ulong(av, o);
    2482             : }
    2483             : 
    2484             : /*Find the exact order of a assuming a^o==1*/
    2485             : GEN
    2486       10656 : Fp_order(GEN a, GEN o, GEN p) {
    2487       10656 :   if (lgefint(p) == 3 && (!o || typ(o) == t_INT))
    2488             :   {
    2489          21 :     ulong pp = p[2], oo = (o && lgefint(o)==3)? uel(o,2): pp-1;
    2490          21 :     return utoi( Fl_order(umodiu(a, pp), oo, pp) );
    2491             :   }
    2492       10635 :   return gen_order(a, o, (void*)p, &Fp_star);
    2493             : }
    2494             : GEN
    2495          70 : Fp_factored_order(GEN a, GEN o, GEN p)
    2496          70 : { return gen_factored_order(a, o, (void*)p, &Fp_star); }
    2497             : 
    2498             : /* return order of a mod p^e, e > 0, pe = p^e */
    2499             : static GEN
    2500          70 : Zp_order(GEN a, GEN p, long e, GEN pe)
    2501             : {
    2502             :   GEN ap, op;
    2503          70 :   if (absequaliu(p, 2))
    2504             :   {
    2505          56 :     if (e == 1) return gen_1;
    2506          56 :     if (e == 2) return mod4(a) == 1? gen_1: gen_2;
    2507          49 :     if (mod4(a) == 1) op = gen_1; else { op = gen_2; a = Fp_sqr(a, pe); }
    2508             :   } else {
    2509          14 :     ap = (e == 1)? a: remii(a,p);
    2510          14 :     op = Fp_order(ap, subiu(p,1), p);
    2511          14 :     if (e == 1) return op;
    2512           0 :     a = Fp_pow(a, op, pe); /* 1 mod p */
    2513             :   }
    2514          49 :   if (equali1(a)) return op;
    2515           7 :   return mulii(op, powiu(p, e - Z_pval(subiu(a,1), p)));
    2516             : }
    2517             : 
    2518             : GEN
    2519          63 : znorder(GEN x, GEN o)
    2520             : {
    2521          63 :   pari_sp av = avma;
    2522             :   GEN b, a;
    2523             : 
    2524          63 :   if (typ(x) != t_INTMOD) pari_err_TYPE("znorder [t_INTMOD expected]",x);
    2525          56 :   b = gel(x,1); a = gel(x,2);
    2526          56 :   if (!equali1(gcdii(a,b))) pari_err_COPRIME("znorder", a,b);
    2527          49 :   if (!o)
    2528             :   {
    2529          35 :     GEN fa = Z_factor(b), P = gel(fa,1), E = gel(fa,2);
    2530          35 :     long i, l = lg(P);
    2531          35 :     o = gen_1;
    2532          70 :     for (i = 1; i < l; i++)
    2533             :     {
    2534          35 :       GEN p = gel(P,i);
    2535          35 :       long e = itos(gel(E,i));
    2536             : 
    2537          35 :       if (l == 2)
    2538          35 :         o = Zp_order(a, p, e, b);
    2539             :       else {
    2540           0 :         GEN pe = powiu(p,e);
    2541           0 :         o = lcmii(o, Zp_order(remii(a,pe), p, e, pe));
    2542             :       }
    2543             :     }
    2544          35 :     return gerepileuptoint(av, o);
    2545             :   }
    2546          14 :   return Fp_order(a, o, b);
    2547             : }
    2548             : 
    2549             : /*********************************************************************/
    2550             : /**               DISCRETE LOGARITHM  in  (Z/nZ)*                   **/
    2551             : /*********************************************************************/
    2552             : static GEN
    2553       58305 : Fp_log_halfgcd(ulong bnd, GEN C, GEN g, GEN p)
    2554             : {
    2555       58305 :   pari_sp av = avma;
    2556             :   GEN h1, h2, F, G;
    2557       58305 :   if (!Fp_ratlift(g,p,C,shifti(C,-1),&h1,&h2)) return gc_NULL(av);
    2558       35013 :   if ((F = Z_issmooth_fact(h1, bnd)) && (G = Z_issmooth_fact(h2, bnd)))
    2559             :   {
    2560         126 :     GEN M = cgetg(3, t_MAT);
    2561         126 :     gel(M,1) = vecsmall_concat(gel(F, 1),gel(G, 1));
    2562         126 :     gel(M,2) = vecsmall_concat(gel(F, 2),zv_neg_inplace(gel(G, 2)));
    2563         126 :     return gerepileupto(av, M);
    2564             :   }
    2565       34887 :   return gc_NULL(av);
    2566             : }
    2567             : 
    2568             : static GEN
    2569       58305 : Fp_log_find_rel(GEN b, ulong bnd, GEN C, GEN p, GEN *g, long *e)
    2570             : {
    2571             :   GEN rel;
    2572       58305 :   do { (*e)++; *g = Fp_mul(*g, b, p); rel = Fp_log_halfgcd(bnd, C, *g, p); }
    2573       58305 :   while (!rel);
    2574         126 :   return rel;
    2575             : }
    2576             : 
    2577             : struct Fp_log_rel
    2578             : {
    2579             :   GEN rel;
    2580             :   ulong prmax;
    2581             :   long nbrel, nbmax, nbgen;
    2582             : };
    2583             : 
    2584             : /* add u^e */
    2585             : static void
    2586        3157 : addifsmooth1(struct Fp_log_rel *r, GEN z, long u, long e)
    2587             : {
    2588        3157 :   pari_sp av = avma;
    2589        3157 :   long off = r->prmax+1;
    2590        3157 :   GEN F = cgetg(3, t_MAT);
    2591        3157 :   gel(F,1) = vecsmall_append(gel(z,1), off+u);
    2592        3157 :   gel(F,2) = vecsmall_append(gel(z,2), e);
    2593        3157 :   gel(r->rel,++r->nbrel) = gerepileupto(av, F);
    2594        3157 : }
    2595             : 
    2596             : /* add u^-1 v^-1 */
    2597             : static void
    2598      104083 : addifsmooth2(struct Fp_log_rel *r, GEN z, long u, long v)
    2599             : {
    2600      104083 :   pari_sp av = avma;
    2601      104083 :   long off = r->prmax+1;
    2602      104083 :   GEN P = mkvecsmall2(off+u,off+v), E = mkvecsmall2(-1,-1);
    2603      104083 :   GEN F = cgetg(3, t_MAT);
    2604      104083 :   gel(F,1) = vecsmall_concat(gel(z,1), P);
    2605      104083 :   gel(F,2) = vecsmall_concat(gel(z,2), E);
    2606      104083 :   gel(r->rel,++r->nbrel) = gerepileupto(av, F);
    2607      104083 : }
    2608             : 
    2609             : /*
    2610             : Let p=C^2+c
    2611             : Solve h = (C+x)*(C+a)-p = 0 [mod l]
    2612             : h= -c+x*(C+a)+C*a = 0  [mod l]
    2613             : x = (c-C*a)/(C+a) [mod l]
    2614             : h = -c+C*(x+a)+a*x
    2615             : */
    2616             : 
    2617             : GEN
    2618       40831 : Fp_log_sieve_worker(long a, long prmax, GEN C, GEN c, GEN Ci, GEN ci, GEN pi, GEN sz)
    2619             : {
    2620       40831 :   pari_sp ltop = avma;
    2621       40831 :   long i, j, th, n = lg(pi)-1, rel = 1;
    2622       40831 :   GEN sieve = zero_zv(a+2)+1;
    2623       40866 :   GEN L = cgetg(1+a+2, t_VEC);
    2624       40853 :   pari_sp av = avma;
    2625       40853 :   GEN z, h = addis(C,a);
    2626       40818 :   if ((z = Z_issmooth_fact(h, prmax)))
    2627             :   {
    2628        3006 :     gel(L, rel++) = mkvec2(z, mkvecsmall3(1, a, -1));
    2629        3014 :     av = avma;
    2630             :   }
    2631    16818295 :   for (i=1; i<=n; i++)
    2632             :   {
    2633    16802230 :     ulong li = pi[i], s = sz[i], al = a % li;
    2634    16802230 :     ulong u, iv = Fl_invsafe(Fl_add(Ci[i],al,li),li);
    2635    17102966 :     if (!iv) continue;
    2636    16679079 :     u = Fl_mul(Fl_sub(ci[i],Fl_mul(Ci[i],al,li),li), iv ,li);
    2637    77374746 :     for(j = u; j<=a; j+=li) sieve[j] += s;
    2638             :   }
    2639       35046 :   if (a)
    2640             :   {
    2641       40740 :     long e = expi(mulis(C,a));
    2642       40777 :     th = e - expu(e) - 1;
    2643          54 :   } else th = -1;
    2644    27986456 :   for (j=0; j<a; j++)
    2645    27944026 :     if (sieve[j]>=th)
    2646             :     {
    2647      119437 :       GEN h = addiu(subii(muliu(C,a+j),c), a*j);
    2648      119343 :       if ((z = Z_issmooth_fact(h, prmax)))
    2649             :       {
    2650      109745 :         gel(L, rel++) = mkvec2(z, mkvecsmall3(2, a, j));
    2651      110003 :         av = avma;
    2652        9423 :       } else set_avma(av);
    2653             :     }
    2654             :   /* j = a */
    2655       42430 :   if (sieve[a]>=th)
    2656             :   {
    2657         476 :     GEN h = addiu(subii(muliu(C,2*a),c), a*a);
    2658         476 :     if ((z = Z_issmooth_fact(h, prmax)))
    2659         385 :       gel(L, rel++) = mkvec2(z, mkvecsmall3(1, a, -2));
    2660             :   }
    2661       42430 :   setlg(L, rel); return gerepilecopy(ltop, L);
    2662             : }
    2663             : 
    2664             : static long
    2665          63 : Fp_log_sieve(struct Fp_log_rel *r, GEN C, GEN c, GEN Ci, GEN ci, GEN pi, GEN sz)
    2666             : {
    2667             :   struct pari_mt pt;
    2668          63 :   long i, j, nb = 0;
    2669          63 :   GEN worker = snm_closure(is_entry("_Fp_log_sieve_worker"),
    2670             :                mkvecn(7, utoi(r->prmax), C, c, Ci, ci, pi, sz));
    2671          63 :   long running, pending = 0;
    2672          63 :   GEN W = zerovec(r->nbgen);
    2673          63 :   mt_queue_start_lim(&pt, worker, r->nbgen);
    2674       41229 :   for (i = 0; (running = (i < r->nbgen)) || pending; i++)
    2675             :   {
    2676             :     GEN done;
    2677             :     long idx;
    2678       41166 :     mt_queue_submit(&pt, i, running ? mkvec(stoi(i)): NULL);
    2679       41166 :     done = mt_queue_get(&pt, &idx, &pending);
    2680       41166 :     if (!done || lg(done)==1) continue;
    2681       35917 :     gel(W, idx+1) = done;
    2682       35917 :     nb += lg(done)-1;
    2683       35917 :     if (DEBUGLEVEL && (i&127)==0)
    2684           0 :       err_printf("%ld%% ",100*nb/r->nbmax);
    2685             :   }
    2686          63 :   mt_queue_end(&pt);
    2687       39550 :   for(j = 1; j <= r->nbgen && r->nbrel < r->nbmax; j++)
    2688             :   {
    2689             :     long ll, m;
    2690       39487 :     GEN L = gel(W,j);
    2691       39487 :     if (isintzero(L)) continue;
    2692       34531 :     ll = lg(L);
    2693      141771 :     for (m=1; m<ll && r->nbrel < r->nbmax ; m++)
    2694             :     {
    2695      107240 :       GEN Lm = gel(L,m), h = gel(Lm, 1), v = gel(Lm, 2);
    2696      107240 :       if (v[1] == 1)
    2697        3157 :         addifsmooth1(r, h, v[2], v[3]);
    2698             :       else
    2699      104083 :         addifsmooth2(r, h, v[2], v[3]);
    2700             :     }
    2701             :   }
    2702          63 :   return j;
    2703             : }
    2704             : 
    2705             : static GEN
    2706         735 : ECP_psi(GEN x, GEN y)
    2707             : {
    2708         735 :   long prec = realprec(x);
    2709         735 :   GEN lx = glog(x, prec), ly = glog(y, prec);
    2710         735 :   GEN u = gdiv(lx, ly);
    2711         735 :   return gpow(u, gneg(u),prec);
    2712             : }
    2713             : 
    2714             : struct computeG
    2715             : {
    2716             :   GEN C;
    2717             :   long bnd, nbi;
    2718             : };
    2719             : 
    2720             : static GEN
    2721         735 : _computeG(void *E, GEN gen)
    2722             : {
    2723         735 :   struct computeG * d = (struct computeG *) E;
    2724         735 :   GEN ps = ECP_psi(gmul(gen,d->C), stoi(d->bnd));
    2725         735 :   return gsub(gmul(gsqr(gen),ps),gmul2n(gaddgs(gen,d->nbi),2));
    2726             : }
    2727             : 
    2728             : static long
    2729          63 : compute_nbgen(GEN C, long bnd, long nbi)
    2730             : {
    2731             :   struct computeG d;
    2732          63 :   d.C = shifti(C, 1);
    2733          63 :   d.bnd = bnd;
    2734          63 :   d.nbi = nbi;
    2735          63 :   return itos(ground(zbrent((void*)&d, _computeG, gen_2, stoi(bnd), DEFAULTPREC)));
    2736             : }
    2737             : 
    2738             : static GEN
    2739        1714 : _psi(void*E, GEN y)
    2740             : {
    2741        1714 :   GEN lx = (GEN) E;
    2742        1714 :   long prec = realprec(lx);
    2743        1714 :   GEN ly = glog(y, prec);
    2744        1714 :   GEN u = gdiv(lx, ly);
    2745        1714 :   return gsub(gdiv(y ,ly), gpow(u, u, prec));
    2746             : }
    2747             : 
    2748             : static GEN
    2749          63 : opt_param(GEN x, long prec)
    2750             : {
    2751          63 :   return zbrent((void*)glog(x,prec), _psi, gen_2, x, prec);
    2752             : }
    2753             : 
    2754             : static GEN
    2755          63 : check_kernel(long nbg, long N, long prmax, GEN C, GEN M, GEN p, GEN m)
    2756             : {
    2757          63 :   pari_sp av = avma;
    2758          63 :   long lM = lg(M)-1, nbcol = lM;
    2759          63 :   long tbs = maxss(1, expu(nbg/expi(m)));
    2760             :   for (;;)
    2761          14 :   {
    2762          77 :     GEN K = FpMs_leftkernel_elt_col(M, nbcol, N, m);
    2763             :     GEN tab;
    2764          77 :     long i, f=0;
    2765          77 :     long l = lg(K), lm = lgefint(m);
    2766          77 :     GEN idx = diviiexact(subiu(p,1),m), g;
    2767             :     pari_timer ti;
    2768          77 :     if (DEBUGLEVEL) timer_start(&ti);
    2769         154 :     for(i=1; i<l; i++)
    2770         154 :       if (signe(gel(K,i)))
    2771          77 :         break;
    2772          77 :     g = Fp_pow(utoi(i), idx, p);
    2773          77 :     tab = Fp_pow_init(g, p, tbs, p);
    2774          77 :     K = FpC_Fp_mul(K, Fp_inv(gel(K,i), m), m);
    2775      128464 :     for(i=1; i<l; i++)
    2776             :     {
    2777      128387 :       GEN k = gel(K,i);
    2778      128387 :       GEN j = i<=prmax ? utoi(i): addis(C,i-(prmax+1));
    2779      128387 :       if (signe(k)==0 || !equalii(Fp_pow_table(tab, k, p), Fp_pow(j, idx, p)))
    2780       76391 :         gel(K,i) = cgetineg(lm);
    2781             :       else
    2782       51996 :         f++;
    2783             :     }
    2784          77 :     if (DEBUGLEVEL) timer_printf(&ti,"found %ld/%ld logs", f, nbg);
    2785          77 :     if(f > (nbg>>1)) return gerepileupto(av, K);
    2786        4585 :     for(i=1; i<=nbcol; i++)
    2787             :     {
    2788        4571 :       long a = 1+random_Fl(lM);
    2789        4571 :       swap(gel(M,a),gel(M,i));
    2790             :     }
    2791          14 :     if (4*nbcol>5*nbg) nbcol = nbcol*9/10;
    2792             :   }
    2793             : }
    2794             : 
    2795             : static GEN
    2796         126 : Fp_log_find_ind(GEN a, GEN K, long prmax, GEN C, GEN p, GEN m)
    2797             : {
    2798         126 :   pari_sp av=avma;
    2799         126 :   GEN aa = gen_1;
    2800         126 :   long AV = 0;
    2801             :   for(;;)
    2802           0 :   {
    2803         126 :     GEN A = Fp_log_find_rel(a, prmax, C, p, &aa, &AV);
    2804         126 :     GEN F = gel(A,1), E = gel(A,2);
    2805         126 :     GEN Ao = gen_0;
    2806         126 :     long i, l = lg(F);
    2807         953 :     for(i=1; i<l; i++)
    2808             :     {
    2809         827 :       GEN Ki = gel(K,F[i]);
    2810         827 :       if (signe(Ki)<0) break;
    2811         827 :       Ao = addii(Ao, mulis(Ki, E[i]));
    2812             :     }
    2813         126 :     if (i==l) return Fp_divu(Ao, AV, m);
    2814           0 :     aa = gerepileuptoint(av, aa);
    2815             :   }
    2816             : }
    2817             : 
    2818             : static GEN
    2819          63 : Fp_log_index(GEN a, GEN b, GEN m, GEN p)
    2820             : {
    2821          63 :   pari_sp av = avma, av2;
    2822          63 :   long i, j, nbi, nbr = 0, nbrow, nbg;
    2823             :   GEN C, c, Ci, ci, pi, pr, sz, l, Ao, Bo, K, d, p_1;
    2824             :   pari_timer ti;
    2825             :   struct Fp_log_rel r;
    2826          63 :   ulong bnds = itou(roundr_safe(opt_param(sqrti(p),DEFAULTPREC)));
    2827          63 :   ulong bnd = 4*bnds;
    2828          63 :   if (!bnds || cmpii(sqru(bnds),m)>=0) return NULL;
    2829             : 
    2830          63 :   p_1 = subiu(p,1);
    2831          63 :   if (!is_pm1(gcdii(m,diviiexact(p_1,m))))
    2832           0 :     m = diviiexact(p_1, Z_ppo(p_1, m));
    2833          63 :   pr = primes_upto_zv(bnd);
    2834          63 :   nbi = lg(pr)-1;
    2835          63 :   C = sqrtremi(p, &c);
    2836          63 :   av2 = avma;
    2837       12796 :   for (i = 1; i <= nbi; ++i)
    2838             :   {
    2839       12733 :     ulong lp = pr[i];
    2840       26894 :     while (lp <= bnd)
    2841             :     {
    2842       14161 :       nbr++;
    2843       14161 :       lp *= pr[i];
    2844             :     }
    2845             :   }
    2846          63 :   pi = cgetg(nbr+1,t_VECSMALL);
    2847          63 :   Ci = cgetg(nbr+1,t_VECSMALL);
    2848          63 :   ci = cgetg(nbr+1,t_VECSMALL);
    2849          63 :   sz = cgetg(nbr+1,t_VECSMALL);
    2850       12796 :   for (i = 1, j = 1; i <= nbi; ++i)
    2851             :   {
    2852       12733 :     ulong lp = pr[i], sp = expu(2*lp-1);
    2853       26894 :     while (lp <= bnd)
    2854             :     {
    2855       14161 :       pi[j] = lp;
    2856       14161 :       Ci[j] = umodiu(C, lp);
    2857       14161 :       ci[j] = umodiu(c, lp);
    2858       14161 :       sz[j] = sp;
    2859       14161 :       lp *= pr[i];
    2860       14161 :       j++;
    2861             :     }
    2862             :   }
    2863          63 :   r.nbrel = 0;
    2864          63 :   r.nbgen = compute_nbgen(C, bnd, nbi);
    2865          63 :   r.nbmax = 2*(nbi+r.nbgen);
    2866          63 :   r.rel = cgetg(r.nbmax+1,t_VEC);
    2867          63 :   r.prmax = pr[nbi];
    2868          63 :   if (DEBUGLEVEL)
    2869             :   {
    2870           0 :     err_printf("bnd=%lu Size FB=%ld extra gen=%ld \n", bnd, nbi, r.nbgen);
    2871           0 :     timer_start(&ti);
    2872             :   }
    2873          63 :   nbg = Fp_log_sieve(&r, C, c, Ci, ci, pi, sz);
    2874          63 :   nbrow = r.prmax + nbg;
    2875          63 :   if (DEBUGLEVEL)
    2876             :   {
    2877           0 :     err_printf("\n");
    2878           0 :     timer_printf(&ti," %ld relations, %ld generators", r.nbrel, nbi+nbg);
    2879             :   }
    2880          63 :   setlg(r.rel,r.nbrel+1);
    2881          63 :   r.rel = gerepilecopy(av2, r.rel);
    2882          63 :   K = check_kernel(nbi+nbrow-r.prmax, nbrow, r.prmax, C, r.rel, p, m);
    2883          63 :   if (DEBUGLEVEL) timer_start(&ti);
    2884          63 :   Ao = Fp_log_find_ind(a, K, r.prmax, C, p, m);
    2885          63 :   if (DEBUGLEVEL) timer_printf(&ti," log element");
    2886          63 :   Bo = Fp_log_find_ind(b, K, r.prmax, C, p, m);
    2887          63 :   if (DEBUGLEVEL) timer_printf(&ti," log generator");
    2888          63 :   d = gcdii(Ao,Bo);
    2889          63 :   l = Fp_div(diviiexact(Ao, d) ,diviiexact(Bo, d), m);
    2890          63 :   if (!equalii(a,Fp_pow(b,l,p))) pari_err_BUG("Fp_log_index");
    2891          63 :   return gerepileuptoint(av, l);
    2892             : }
    2893             : 
    2894             : static int
    2895     4642731 : Fp_log_use_index(long e, long p)
    2896             : {
    2897     4642731 :   return (e >= 27 && 20*(p+6)<=e*e);
    2898             : }
    2899             : 
    2900             : /* Trivial cases a = 1, -1. Return x s.t. g^x = a or [] if no such x exist */
    2901             : static GEN
    2902     8442668 : Fp_easylog(void *E, GEN a, GEN g, GEN ord)
    2903             : {
    2904     8442668 :   pari_sp av = avma;
    2905     8442668 :   GEN p = (GEN)E;
    2906             :   /* assume a reduced mod p, p not necessarily prime */
    2907     8442668 :   if (equali1(a)) return gen_0;
    2908             :   /* p > 2 */
    2909     5421247 :   if (equalii(subiu(p,1), a))  /* -1 */
    2910             :   {
    2911             :     pari_sp av2;
    2912             :     GEN t;
    2913     1321379 :     ord = get_arith_Z(ord);
    2914     1321379 :     if (mpodd(ord)) { set_avma(av); return cgetg(1, t_VEC); } /* no solution */
    2915     1321365 :     t = shifti(ord,-1); /* only possible solution */
    2916     1321364 :     av2 = avma;
    2917     1321364 :     if (!equalii(Fp_pow(g, t, p), a)) { set_avma(av); return cgetg(1, t_VEC); }
    2918     1321335 :     set_avma(av2); return gerepileuptoint(av, t);
    2919             :   }
    2920     4099870 :   if (typ(ord)==t_INT && BPSW_psp(p) && Fp_log_use_index(expi(ord),expi(p)))
    2921          63 :     return Fp_log_index(a, g, ord, p);
    2922     4099807 :   return gc_NULL(av); /* not easy */
    2923             : }
    2924             : 
    2925             : GEN
    2926     3913070 : Fp_log(GEN a, GEN g, GEN ord, GEN p)
    2927             : {
    2928     3913070 :   GEN v = get_arith_ZZM(ord);
    2929     3913048 :   GEN F = gmael(v,2,1);
    2930     3913048 :   long lF = lg(F)-1, lmax;
    2931     3913048 :   if (lF == 0) return equali1(a)? gen_0: cgetg(1, t_VEC);
    2932     3913020 :   lmax = expi(gel(F,lF));
    2933     3913022 :   if (BPSW_psp(p) && Fp_log_use_index(lmax,expi(p)))
    2934          91 :     v = mkvec2(gel(v,1),ZM_famat_limit(gel(v,2),int2n(27)));
    2935     3913018 :   return gen_PH_log(a,g,v,(void*)p,&Fp_star);
    2936             : }
    2937             : 
    2938             : /* assume !(p & HIGHMASK) */
    2939             : static ulong
    2940      130283 : Fl_log_naive(ulong a, ulong g, ulong ord, ulong p)
    2941             : {
    2942      130283 :   ulong i, h=1;
    2943      357162 :   for (i = 0; i < ord; i++, h = (h * g) % p)
    2944      357162 :     if (a==h) return i;
    2945           0 :   return ~0UL;
    2946             : }
    2947             : 
    2948             : static ulong
    2949       22215 : Fl_log_naive_pre(ulong a, ulong g, ulong ord, ulong p, ulong pi)
    2950             : {
    2951       22215 :   ulong i, h=1;
    2952       56012 :   for (i = 0; i < ord; i++, h = Fl_mul_pre(h, g, p, pi))
    2953       56012 :     if (a==h) return i;
    2954           0 :   return ~0UL;
    2955             : }
    2956             : 
    2957             : static ulong
    2958           0 : Fl_log_Fp(ulong a, ulong g, ulong ord, ulong p)
    2959             : {
    2960           0 :   pari_sp av = avma;
    2961           0 :   GEN r = Fp_log(utoi(a),utoi(g),utoi(ord),utoi(p));
    2962           0 :   return gc_ulong(av, typ(r)==t_INT ? itou(r): ~0UL);
    2963             : }
    2964             : 
    2965             : /* allow pi = 0 */
    2966             : ulong
    2967       22593 : Fl_log_pre(ulong a, ulong g, ulong ord, ulong p, ulong pi)
    2968             : {
    2969       22593 :   if (!pi) return Fl_log(a, g, ord, p);
    2970       22215 :   if (ord <= 200) return Fl_log_naive_pre(a, g, ord, p, pi);
    2971           0 :   return Fl_log_Fp(a, g, ord, p);
    2972             : }
    2973             : 
    2974             : ulong
    2975      130283 : Fl_log(ulong a, ulong g, ulong ord, ulong p)
    2976             : {
    2977      130283 :   if (ord <= 200)
    2978           0 :     return (p&HIGHMASK)? Fl_log_naive_pre(a, g, ord, p, get_Fl_red(p))
    2979      130283 :                        : Fl_log_naive(a, g, ord, p);
    2980           0 :   return Fl_log_Fp(a, g, ord, p);
    2981             : }
    2982             : 
    2983             : /* find x such that h = g^x mod N > 1, N = prod_{i <= l} P[i]^E[i], P[i] prime.
    2984             :  * PHI[l] = eulerphi(N / P[l]^E[l]).   Destroys P/E */
    2985             : static GEN
    2986         126 : znlog_rec(GEN h, GEN g, GEN N, GEN P, GEN E, GEN PHI)
    2987             : {
    2988         126 :   long l = lg(P) - 1, e = E[l];
    2989         126 :   GEN p = gel(P, l), phi = gel(PHI,l), pe = e == 1? p: powiu(p, e);
    2990             :   GEN a,b, hp,gp, hpe,gpe, ogpe; /* = order(g mod p^e) | p^(e-1)(p-1) */
    2991             : 
    2992         126 :   if (l == 1) {
    2993          98 :     hpe = h;
    2994          98 :     gpe = g;
    2995             :   } else {
    2996          28 :     hpe = modii(h, pe);
    2997          28 :     gpe = modii(g, pe);
    2998             :   }
    2999         126 :   if (e == 1) {
    3000          42 :     hp = hpe;
    3001          42 :     gp = gpe;
    3002             :   } else {
    3003          84 :     hp = remii(hpe, p);
    3004          84 :     gp = remii(gpe, p);
    3005             :   }
    3006         126 :   if (hp == gen_0 || gp == gen_0) return NULL;
    3007         105 :   if (absequaliu(p, 2))
    3008             :   {
    3009          35 :     GEN n = int2n(e);
    3010          35 :     ogpe = Zp_order(gpe, gen_2, e, n);
    3011          35 :     a = Fp_log(hpe, gpe, ogpe, n);
    3012          35 :     if (typ(a) != t_INT) return NULL;
    3013             :   }
    3014             :   else
    3015             :   { /* Avoid black box groups: (Z/p^2)^* / (Z/p)^* ~ (Z/pZ, +), where DL
    3016             :        is trivial */
    3017             :     /* [order(gp), factor(order(gp))] */
    3018          70 :     GEN v = Fp_factored_order(gp, subiu(p,1), p);
    3019          70 :     GEN ogp = gel(v,1);
    3020          70 :     if (!equali1(Fp_pow(hp, ogp, p))) return NULL;
    3021          70 :     a = Fp_log(hp, gp, v, p);
    3022          70 :     if (typ(a) != t_INT) return NULL;
    3023          70 :     if (e == 1) ogpe = ogp;
    3024             :     else
    3025             :     { /* find a s.t. g^a = h (mod p^e), p odd prime, e > 0, (h,p) = 1 */
    3026             :       /* use p-adic log: O(log p + e) mul*/
    3027             :       long vpogpe, vpohpe;
    3028             : 
    3029          28 :       hpe = Fp_mul(hpe, Fp_pow(gpe, negi(a), pe), pe);
    3030          28 :       gpe = Fp_pow(gpe, ogp, pe);
    3031             :       /* g,h = 1 mod p; compute b s.t. h = g^b */
    3032             : 
    3033             :       /* v_p(order g mod pe) */
    3034          28 :       vpogpe = equali1(gpe)? 0: e - Z_pval(subiu(gpe,1), p);
    3035             :       /* v_p(order h mod pe) */
    3036          28 :       vpohpe = equali1(hpe)? 0: e - Z_pval(subiu(hpe,1), p);
    3037          28 :       if (vpohpe > vpogpe) return NULL;
    3038             : 
    3039          28 :       ogpe = mulii(ogp, powiu(p, vpogpe)); /* order g mod p^e */
    3040          28 :       if (is_pm1(gpe)) return is_pm1(hpe)? a: NULL;
    3041          28 :       b = gdiv(Qp_log(cvtop(hpe, p, e)), Qp_log(cvtop(gpe, p, e)));
    3042          28 :       a = addii(a, mulii(ogp, padic_to_Q(b)));
    3043             :     }
    3044             :   }
    3045             :   /* gp^a = hp => x = a mod ogpe => generalized Pohlig-Hellman strategy */
    3046          91 :   if (l == 1) return a;
    3047             : 
    3048          28 :   N = diviiexact(N, pe); /* make N coprime to p */
    3049          28 :   h = Fp_mul(h, Fp_pow(g, modii(negi(a), phi), N), N);
    3050          28 :   g = Fp_pow(g, modii(ogpe, phi), N);
    3051          28 :   setlg(P, l); /* remove last element */
    3052          28 :   setlg(E, l);
    3053          28 :   b = znlog_rec(h, g, N, P, E, PHI);
    3054          28 :   if (!b) return NULL;
    3055          28 :   return addmulii(a, b, ogpe);
    3056             : }
    3057             : 
    3058             : static GEN
    3059          98 : get_PHI(GEN P, GEN E)
    3060             : {
    3061          98 :   long i, l = lg(P);
    3062          98 :   GEN PHI = cgetg(l, t_VEC);
    3063          98 :   gel(PHI,1) = gen_1;
    3064         126 :   for (i=1; i<l-1; i++)
    3065             :   {
    3066          28 :     GEN t, p = gel(P,i);
    3067          28 :     long e = E[i];
    3068          28 :     t = mulii(powiu(p, e-1), subiu(p,1));
    3069          28 :     if (i > 1) t = mulii(t, gel(PHI,i));
    3070          28 :     gel(PHI,i+1) = t;
    3071             :   }
    3072          98 :   return PHI;
    3073             : }
    3074             : 
    3075             : GEN
    3076         238 : znlog(GEN h, GEN g, GEN o)
    3077             : {
    3078         238 :   pari_sp av = avma;
    3079             :   GEN N, fa, P, E, x;
    3080         238 :   switch (typ(g))
    3081             :   {
    3082          28 :     case t_PADIC:
    3083             :     {
    3084          28 :       GEN p = gel(g,2);
    3085          28 :       long v = valp(g);
    3086          28 :       if (v < 0) pari_err_DIM("znlog");
    3087          28 :       if (v > 0) {
    3088           0 :         long k = gvaluation(h, p);
    3089           0 :         if (k % v) return cgetg(1,t_VEC);
    3090           0 :         k /= v;
    3091           0 :         if (!gequal(h, gpowgs(g,k))) { set_avma(av); return cgetg(1,t_VEC); }
    3092           0 :         return gc_stoi(av, k);
    3093             :       }
    3094          28 :       N = gel(g,3);
    3095          28 :       g = Rg_to_Fp(g, N);
    3096          28 :       break;
    3097             :     }
    3098         203 :     case t_INTMOD:
    3099         203 :       N = gel(g,1);
    3100         203 :       g = gel(g,2); break;
    3101           7 :     default: pari_err_TYPE("znlog", g);
    3102             :       return NULL; /* LCOV_EXCL_LINE */
    3103             :   }
    3104         231 :   if (equali1(N)) { set_avma(av); return gen_0; }
    3105         231 :   h = Rg_to_Fp(h, N);
    3106         224 :   if (o) return gerepileupto(av, Fp_log(h, g, o, N));
    3107          98 :   fa = Z_factor(N);
    3108          98 :   P = gel(fa,1);
    3109          98 :   E = vec_to_vecsmall(gel(fa,2));
    3110          98 :   x = znlog_rec(h, g, N, P, E, get_PHI(P,E));
    3111          98 :   if (!x) { set_avma(av); return cgetg(1,t_VEC); }
    3112          63 :   return gerepileuptoint(av, x);
    3113             : }
    3114             : 
    3115             : GEN
    3116       62832 : Fp_sqrtn(GEN a, GEN n, GEN p, GEN *zeta)
    3117             : {
    3118       62832 :   if (lgefint(p)==3)
    3119             :   {
    3120       62398 :     long nn = itos_or_0(n);
    3121       62398 :     if (nn)
    3122             :     {
    3123       62398 :       ulong pp = p[2];
    3124             :       ulong uz;
    3125       62398 :       ulong r = Fl_sqrtn(umodiu(a,pp),nn,pp, zeta ? &uz:NULL);
    3126       62377 :       if (r==ULONG_MAX) return NULL;
    3127       62321 :       if (zeta) *zeta = utoi(uz);
    3128       62321 :       return utoi(r);
    3129             :     }
    3130             :   }
    3131         434 :   a = modii(a,p);
    3132         434 :   if (!signe(a))
    3133             :   {
    3134           0 :     if (zeta) *zeta = gen_1;
    3135           0 :     if (signe(n) < 0) pari_err_INV("Fp_sqrtn", mkintmod(gen_0,p));
    3136           0 :     return gen_0;
    3137             :   }
    3138         434 :   if (absequaliu(n,2))
    3139             :   {
    3140         238 :     if (zeta) *zeta = subiu(p,1);
    3141         238 :     return signe(n) > 0 ? Fp_sqrt(a,p): Fp_sqrt(Fp_inv(a, p),p);
    3142             :   }
    3143         196 :   return gen_Shanks_sqrtn(a,n,subiu(p,1),zeta,(void*)p,&Fp_star);
    3144             : }
    3145             : 
    3146             : /*********************************************************************/
    3147             : /**                              FACTORIAL                          **/
    3148             : /*********************************************************************/
    3149             : GEN
    3150       84157 : mulu_interval_step(ulong a, ulong b, ulong step)
    3151             : {
    3152       84157 :   pari_sp av = avma;
    3153             :   ulong k, l, N, n;
    3154             :   long lx;
    3155             :   GEN x;
    3156             : 
    3157       84157 :   if (!a) return gen_0;
    3158       84157 :   if (step == 1) return mulu_interval(a, b);
    3159       84157 :   n = 1 + (b-a) / step;
    3160       84157 :   b -= (b-a) % step;
    3161       84157 :   if (n < 61)
    3162             :   {
    3163       82708 :     if (n == 1) return utoipos(a);
    3164       63539 :     x = muluu(a,a+step); if (n == 2) return x;
    3165      511199 :     for (k=a+2*step; k<=b; k+=step) x = mului(k,x);
    3166       49656 :     return gerepileuptoint(av, x);
    3167             :   }
    3168             :   /* step | b-a */
    3169        1449 :   lx = 1; x = cgetg(2 + n/2, t_VEC);
    3170        1449 :   N = b + a;
    3171        1449 :   for (k = a;; k += step)
    3172             :   {
    3173      230241 :     l = N - k; if (l <= k) break;
    3174      228792 :     gel(x,lx++) = muluu(k,l);
    3175             :   }
    3176        1449 :   if (l == k) gel(x,lx++) = utoipos(k);
    3177        1449 :   setlg(x, lx);
    3178        1449 :   return gerepileuptoint(av, ZV_prod(x));
    3179             : }
    3180             : /* return a * (a+1) * ... * b. Assume a <= b  [ note: factoring out powers of 2
    3181             :  * first is slower ... ] */
    3182             : GEN
    3183      156497 : mulu_interval(ulong a, ulong b)
    3184             : {
    3185      156497 :   pari_sp av = avma;
    3186             :   ulong k, l, N, n;
    3187             :   long lx;
    3188             :   GEN x;
    3189             : 
    3190      156497 :   if (!a) return gen_0;
    3191      156497 :   n = b - a + 1;
    3192      156497 :   if (n < 61)
    3193             :   {
    3194      156497 :     if (n == 1) return utoipos(a);
    3195      106174 :     x = muluu(a,a+1); if (n == 2) return x;
    3196      371718 :     for (k=a+2; k<b; k++) x = mului(k,x);
    3197             :     /* avoid k <= b: broken if b = ULONG_MAX */
    3198       92132 :     return gerepileuptoint(av, mului(b,x));
    3199             :   }
    3200           0 :   lx = 1; x = cgetg(2 + n/2, t_VEC);
    3201           0 :   N = b + a;
    3202           0 :   for (k = a;; k++)
    3203             :   {
    3204           0 :     l = N - k; if (l <= k) break;
    3205           0 :     gel(x,lx++) = muluu(k,l);
    3206             :   }
    3207           0 :   if (l == k) gel(x,lx++) = utoipos(k);
    3208           0 :   setlg(x, lx);
    3209           0 :   return gerepileuptoint(av, ZV_prod(x));
    3210             : }
    3211             : GEN
    3212         588 : muls_interval(long a, long b)
    3213             : {
    3214         588 :   pari_sp av = avma;
    3215         588 :   long lx, k, l, N, n = b - a + 1;
    3216             :   GEN x;
    3217             : 
    3218         588 :   if (a <= 0 && b >= 0) return gen_0;
    3219         315 :   if (n < 61)
    3220             :   {
    3221         315 :     x = stoi(a);
    3222         553 :     for (k=a+1; k<=b; k++) x = mulsi(k,x);
    3223         315 :     return gerepileuptoint(av, x);
    3224             :   }
    3225           0 :   lx = 1; x = cgetg(2 + n/2, t_VEC);
    3226           0 :   N = b + a;
    3227           0 :   for (k = a;; k++)
    3228             :   {
    3229           0 :     l = N - k; if (l <= k) break;
    3230           0 :     gel(x,lx++) = mulss(k,l);
    3231             :   }
    3232           0 :   if (l == k) gel(x,lx++) = stoi(k);
    3233           0 :   setlg(x, lx);
    3234           0 :   return gerepileuptoint(av, ZV_prod(x));
    3235             : }
    3236             : 
    3237             : GEN
    3238      525376 : mpfact(long n)
    3239             : {
    3240      525376 :   pari_sp av = avma;
    3241             :   GEN a, v;
    3242             :   long k;
    3243      525376 :   if (n <= 12) switch(n)
    3244             :   {
    3245      446373 :     case 0: case 1: return gen_1;
    3246       36916 :     case 2: return gen_2;
    3247        3390 :     case 3: return utoipos(6);
    3248        4233 :     case 4: return utoipos(24);
    3249        2859 :     case 5: return utoipos(120);
    3250        2516 :     case 6: return utoipos(720);
    3251        2410 :     case 7: return utoipos(5040);
    3252        2411 :     case 8: return utoipos(40320);
    3253        2432 :     case 9: return utoipos(362880);
    3254        2684 :     case 10:return utoipos(3628800);
    3255        1425 :     case 11:return utoipos(39916800);
    3256         579 :     case 12:return utoipos(479001600);
    3257           0 :     default: pari_err_DOMAIN("factorial", "argument","<",gen_0,stoi(n));
    3258             :   }
    3259       17148 :   v = cgetg(expu(n) + 2, t_VEC);
    3260       17148 :   for (k = 1;; k++)
    3261       80363 :   {
    3262       97511 :     long m = n >> (k-1), l;
    3263       97511 :     if (m <= 2) break;
    3264       80363 :     l = (1 + (n >> k)) | 1;
    3265             :     /* product of odd numbers in ]n / 2^k, 2 / 2^(k-1)] */
    3266       80363 :     a = mulu_interval_step(l, m, 2);
    3267       80363 :     gel(v,k) = k == 1? a: powiu(a, k);
    3268             :   }
    3269       80363 :   a = gel(v,--k); while (--k) a = mulii(a, gel(v,k));
    3270       17148 :   a = shifti(a, factorial_lval(n, 2));
    3271       17148 :   return gerepileuptoint(av, a);
    3272             : }
    3273             : 
    3274             : ulong
    3275       56752 : factorial_Fl(long n, ulong p)
    3276             : {
    3277             :   long k;
    3278             :   ulong v;
    3279       56752 :   if (p <= (ulong)n) return 0;
    3280       56752 :   v = Fl_powu(2, factorial_lval(n, 2), p);
    3281       56796 :   for (k = 1;; k++)
    3282      142389 :   {
    3283      199185 :     long m = n >> (k-1), l, i;
    3284      199185 :     ulong a = 1;
    3285      199185 :     if (m <= 2) break;
    3286      142408 :     l = (1 + (n >> k)) | 1;
    3287             :     /* product of odd numbers in ]n / 2^k, 2 / 2^(k-1)] */
    3288      779701 :     for (i=l; i<=m; i+=2)
    3289      637306 :       a = Fl_mul(a, i, p);
    3290      142395 :     v = Fl_mul(v, k == 1? a: Fl_powu(a, k, p), p);
    3291             :   }
    3292       56777 :   return v;
    3293             : }
    3294             : 
    3295             : GEN
    3296         158 : factorial_Fp(long n, GEN p)
    3297             : {
    3298         158 :   pari_sp av = avma;
    3299             :   long k;
    3300         158 :   GEN v = Fp_powu(gen_2, factorial_lval(n, 2), p);
    3301         158 :   for (k = 1;; k++)
    3302         344 :   {
    3303         502 :     long m = n >> (k-1), l, i;
    3304         502 :     GEN a = gen_1;
    3305         502 :     if (m <= 2) break;
    3306         344 :     l = (1 + (n >> k)) | 1;
    3307             :     /* product of odd numbers in ]n / 2^k, 2 / 2^(k-1)] */
    3308         850 :     for (i=l; i<=m; i+=2)
    3309         506 :       a = Fp_mulu(a, i, p);
    3310         344 :     v = Fp_mul(v, k == 1? a: Fp_powu(a, k, p), p);
    3311         344 :     v = gerepileuptoint(av, v);
    3312             :   }
    3313         158 :   return v;
    3314             : }
    3315             : 
    3316             : /*******************************************************************/
    3317             : /**                      LUCAS & FIBONACCI                        **/
    3318             : /*******************************************************************/
    3319             : static void
    3320          56 : lucas(ulong n, GEN *a, GEN *b)
    3321             : {
    3322             :   GEN z, t, zt;
    3323          56 :   if (!n) { *a = gen_2; *b = gen_1; return; }
    3324          49 :   lucas(n >> 1, &z, &t); zt = mulii(z, t);
    3325          49 :   switch(n & 3) {
    3326          14 :     case  0: *a = subiu(sqri(z),2); *b = subiu(zt,1); break;
    3327          14 :     case  1: *a = subiu(zt,1);      *b = addiu(sqri(t),2); break;
    3328           7 :     case  2: *a = addiu(sqri(z),2); *b = addiu(zt,1); break;
    3329          14 :     case  3: *a = addiu(zt,1);      *b = subiu(sqri(t),2);
    3330             :   }
    3331          49 : }
    3332             : 
    3333             : GEN
    3334           7 : fibo(long n)
    3335             : {
    3336           7 :   pari_sp av = avma;
    3337             :   GEN a, b;
    3338           7 :   if (!n) return gen_0;
    3339           7 :   lucas((ulong)(labs(n)-1), &a, &b);
    3340           7 :   a = diviuexact(addii(shifti(a,1),b), 5);
    3341           7 :   if (n < 0 && !odd(n)) setsigne(a, -1);
    3342           7 :   return gerepileuptoint(av, a);
    3343             : }
    3344             : 
    3345             : /*******************************************************************/
    3346             : /*                      CONTINUED FRACTIONS                        */
    3347             : /*******************************************************************/
    3348             : static GEN
    3349     2830748 : icopy_lg(GEN x, long l)
    3350             : {
    3351     2830748 :   long lx = lgefint(x);
    3352             :   GEN y;
    3353             : 
    3354     2830748 :   if (lx >= l) return icopy(x);
    3355          49 :   y = cgeti(l); affii(x, y); return y;
    3356             : }
    3357             : 
    3358             : /* continued fraction of a/b. If y != NULL, stop when partial quotients
    3359             :  * differ from y */
    3360             : static GEN
    3361     2831102 : Qsfcont(GEN a, GEN b, GEN y, ulong k)
    3362             : {
    3363             :   GEN  z, c;
    3364     2831102 :   ulong i, l, ly = lgefint(b);
    3365             : 
    3366             :   /* times 1 / log2( (1+sqrt(5)) / 2 )  */
    3367     2831102 :   l = (ulong)(3 + bit_accuracy_mul(ly, 1.44042009041256));
    3368     2831102 :   if (k > 0 && k+1 > 0 && l > k+1) l = k+1; /* beware overflow */
    3369     2831102 :   if (l > LGBITS) l = LGBITS;
    3370             : 
    3371     2831102 :   z = cgetg(l,t_VEC);
    3372     2831102 :   l--;
    3373     2831102 :   if (y) {
    3374         354 :     pari_sp av = avma;
    3375         354 :     if (l >= (ulong)lg(y)) l = lg(y)-1;
    3376       28713 :     for (i = 1; i <= l; i++)
    3377             :     {
    3378       28482 :       GEN q = gel(y,i);
    3379       28482 :       gel(z,i) = q;
    3380       28482 :       c = b; if (!gequal1(q)) c = mulii(q, b);
    3381       28482 :       c = subii(a, c);
    3382       28482 :       if (signe(c) < 0)
    3383             :       { /* partial quotient too large */
    3384         117 :         c = addii(c, b);
    3385         117 :         if (signe(c) >= 0) i++; /* by 1 */
    3386         117 :         break;
    3387             :       }
    3388       28365 :       if (cmpii(c, b) >= 0)
    3389             :       { /* partial quotient too small */
    3390           6 :         c = subii(c, b);
    3391           6 :         if (cmpii(c, b) < 0) {
    3392             :           /* by 1. If next quotient is 1 in y, add 1 */
    3393           6 :           if (i < l && equali1(gel(y,i+1))) gel(z,i) = addiu(q,1);
    3394           6 :           i++;
    3395             :         }
    3396           6 :         break;
    3397             :       }
    3398       28359 :       if ((i & 0xff) == 0) gerepileall(av, 2, &b, &c);
    3399       28359 :       a = b; b = c;
    3400             :     }
    3401             :   } else {
    3402     2830748 :     a = icopy_lg(a, ly);
    3403     2830748 :     b = icopy(b);
    3404    23451289 :     for (i = 1; i <= l; i++)
    3405             :     {
    3406    23450959 :       gel(z,i) = truedvmdii(a,b,&c);
    3407    23450959 :       if (c == gen_0) { i++; break; }
    3408    20620541 :       affii(c, a); cgiv(c); c = a;
    3409    20620541 :       a = b; b = c;
    3410             :     }
    3411             :   }
    3412     2831102 :   i--;
    3413     2831102 :   if (i > 1 && gequal1(gel(z,i)))
    3414             :   {
    3415         100 :     cgiv(gel(z,i)); --i;
    3416         100 :     gel(z,i) = addui(1, gel(z,i)); /* unclean: leave old z[i] on stack */
    3417             :   }
    3418     2831102 :   setlg(z,i+1); return z;
    3419             : }
    3420             : 
    3421             : static GEN
    3422           0 : sersfcont(GEN a, GEN b, long k)
    3423             : {
    3424           0 :   long i, l = typ(a) == t_POL? lg(a): 3;
    3425             :   GEN y, c;
    3426           0 :   if (lg(b) > l) l = lg(b);
    3427           0 :   if (k > 0 && l > k+1) l = k+1;
    3428           0 :   y = cgetg(l,t_VEC);
    3429           0 :   for (i=1; i<l; i++)
    3430             :   {
    3431           0 :     gel(y,i) = poldivrem(a,b,&c);
    3432           0 :     if (gequal0(c)) { i++; break; }
    3433           0 :     a = b; b = c;
    3434             :   }
    3435           0 :   setlg(y, i); return y;
    3436             : }
    3437             : 
    3438             : GEN
    3439     2831763 : gboundcf(GEN x, long k)
    3440             : {
    3441             :   pari_sp av;
    3442     2831763 :   long tx = typ(x), e;
    3443             :   GEN y, a, b, c;
    3444             : 
    3445     2831763 :   if (k < 0) pari_err_DOMAIN("gboundcf","nmax","<",gen_0,stoi(k));
    3446     2831756 :   if (is_scalar_t(tx))
    3447             :   {
    3448     2831756 :     if (gequal0(x)) return mkvec(gen_0);
    3449     2831651 :     switch(tx)
    3450             :     {
    3451         896 :       case t_INT: return mkveccopy(x);
    3452         361 :       case t_REAL:
    3453         361 :         av = avma;
    3454         361 :         c = mantissa_real(x,&e);
    3455         361 :         if (e < 0) pari_err_PREC("gboundcf");
    3456         354 :         y = int2n(e);
    3457         354 :         a = Qsfcont(c,y, NULL, k);
    3458         354 :         b = addsi(signe(x), c);
    3459         354 :         return gerepilecopy(av, Qsfcont(b,y, a, k));
    3460             : 
    3461     2830394 :       case t_FRAC:
    3462     2830394 :         av = avma;
    3463     2830394 :         return gerepileupto(av, Qsfcont(gel(x,1),gel(x,2), NULL, k));
    3464             :     }
    3465           0 :     pari_err_TYPE("gboundcf",x);
    3466             :   }
    3467             : 
    3468           0 :   switch(tx)
    3469             :   {
    3470           0 :     case t_POL: return mkveccopy(x);
    3471           0 :     case t_SER:
    3472           0 :       av = avma;
    3473           0 :       return gerepileupto(av, gboundcf(ser2rfrac_i(x), k));
    3474           0 :     case t_RFRAC:
    3475           0 :       av = avma;
    3476           0 :       return gerepilecopy(av, sersfcont(gel(x,1), gel(x,2), k));
    3477             :   }
    3478           0 :   pari_err_TYPE("gboundcf",x);
    3479             :   return NULL; /* LCOV_EXCL_LINE */
    3480             : }
    3481             : 
    3482             : static GEN
    3483          14 : sfcont2(GEN b, GEN x, long k)
    3484             : {
    3485          14 :   pari_sp av = avma;
    3486          14 :   long lb = lg(b), tx = typ(x), i;
    3487             :   GEN y,p1;
    3488             : 
    3489          14 :   if (k)
    3490             :   {
    3491           7 :     if (k >= lb) pari_err_DIM("contfrac [too few denominators]");
    3492           0 :     lb = k+1;
    3493             :   }
    3494           7 :   y = cgetg(lb,t_VEC);
    3495           7 :   if (lb==1) return y;
    3496           7 :   if (is_scalar_t(tx))
    3497             :   {
    3498           7 :     if (!is_intreal_t(tx) && tx != t_FRAC) pari_err_TYPE("sfcont2",x);
    3499             :   }
    3500           0 :   else if (tx == t_SER) x = ser2rfrac_i(x);
    3501             : 
    3502           7 :   if (!gequal1(gel(b,1))) x = gmul(gel(b,1),x);
    3503           7 :   for (i = 1;;)
    3504             :   {
    3505          35 :     if (tx == t_REAL)
    3506             :     {
    3507          35 :       long e = expo(x);
    3508          35 :       if (e > 0 && nbits2prec(e+1) > realprec(x)) break;
    3509          35 :       gel(y,i) = floorr(x);
    3510          35 :       p1 = subri(x, gel(y,i));
    3511             :     }
    3512             :     else
    3513             :     {
    3514           0 :       gel(y,i) = gfloor(x);
    3515           0 :       p1 = gsub(x, gel(y,i));
    3516             :     }
    3517          35 :     if (++i >= lb) break;
    3518          28 :     if (gequal0(p1)) break;
    3519          28 :     x = gdiv(gel(b,i),p1);
    3520             :   }
    3521           7 :   setlg(y,i);
    3522           7 :   return gerepilecopy(av,y);
    3523             : }
    3524             : 
    3525             : GEN
    3526         112 : gcf(GEN x) { return gboundcf(x,0); }
    3527             : GEN
    3528           0 : gcf2(GEN b, GEN x) { return contfrac0(x,b,0); }
    3529             : GEN
    3530          49 : contfrac0(GEN x, GEN b, long nmax)
    3531             : {
    3532             :   long tb;
    3533             : 
    3534          49 :   if (!b) return gboundcf(x,nmax);
    3535          28 :   tb = typ(b);
    3536          28 :   if (tb == t_INT) return gboundcf(x,itos(b));
    3537          21 :   if (! is_vec_t(tb)) pari_err_TYPE("contfrac0",b);
    3538          21 :   if (nmax < 0) pari_err_DOMAIN("contfrac","nmax","<",gen_0,stoi(nmax));
    3539          14 :   return sfcont2(b,x,nmax);
    3540             : }
    3541             : 
    3542             : GEN
    3543         252 : contfracpnqn(GEN x, long n)
    3544             : {
    3545         252 :   pari_sp av = avma;
    3546         252 :   long i, lx = lg(x);
    3547             :   GEN M,A,B, p0,p1, q0,q1;
    3548             : 
    3549         252 :   if (lx == 1)
    3550             :   {
    3551          28 :     if (! is_matvec_t(typ(x))) pari_err_TYPE("pnqn",x);
    3552          21 :     if (n >= 0) return cgetg(1,t_MAT);
    3553           7 :     return matid(2);
    3554             :   }
    3555         224 :   switch(typ(x))
    3556             :   {
    3557         182 :     case t_VEC: case t_COL: A = x; B = NULL; break;
    3558          42 :     case t_MAT:
    3559          42 :       switch(lgcols(x))
    3560             :       {
    3561           0 :         case 2: A = row(x,1); B = NULL; break;
    3562          35 :         case 3: A = row(x,2); B = row(x,1); break;
    3563           7 :         default: pari_err_DIM("pnqn [ nbrows != 1,2 ]");
    3564             :                  return NULL; /*LCOV_EXCL_LINE*/
    3565             :       }
    3566          35 :       break;
    3567           0 :     default: pari_err_TYPE("pnqn",x);
    3568             :       return NULL; /*LCOV_EXCL_LINE*/
    3569             :   }
    3570         217 :   p1 = gel(A,1);
    3571         217 :   q1 = B? gel(B,1): gen_1; /* p[0], q[0] */
    3572         217 :   if (n >= 0)
    3573             :   {
    3574         182 :     lx = minss(lx, n+2);
    3575         182 :     if (lx == 2) return gerepilecopy(av, mkmat(mkcol2(p1,q1)));
    3576             :   }
    3577          35 :   else if (lx == 2)
    3578           7 :     return gerepilecopy(av, mkmat2(mkcol2(p1,q1), mkcol2(gen_1,gen_0)));
    3579             :   /* lx >= 3 */
    3580         119 :   p0 = gen_1;
    3581         119 :   q0 = gen_0; /* p[-1], q[-1] */
    3582         119 :   M = cgetg(lx, t_MAT);
    3583         119 :   gel(M,1) = mkcol2(p1,q1);
    3584         399 :   for (i=2; i<lx; i++)
    3585             :   {
    3586         280 :     GEN a = gel(A,i), p2,q2;
    3587         280 :     if (B) {
    3588          84 :       GEN b = gel(B,i);
    3589          84 :       p0 = gmul(b,p0);
    3590          84 :       q0 = gmul(b,q0);
    3591             :     }
    3592         280 :     p2 = gadd(gmul(a,p1),p0); p0=p1; p1=p2;
    3593         280 :     q2 = gadd(gmul(a,q1),q0); q0=q1; q1=q2;
    3594         280 :     gel(M,i) = mkcol2(p1,q1);
    3595             :   }
    3596         119 :   if (n < 0) M = mkmat2(gel(M,lx-1), gel(M,lx-2));
    3597         119 :   return gerepilecopy(av, M);
    3598             : }
    3599             : GEN
    3600           0 : pnqn(GEN x) { return contfracpnqn(x,-1); }
    3601             : /* x = [a0, ..., an] from gboundcf, n >= 0;
    3602             :  * return [[p0, ..., pn], [q0,...,qn]] */
    3603             : GEN
    3604      609308 : ZV_allpnqn(GEN x)
    3605             : {
    3606      609308 :   long i, lx = lg(x);
    3607      609308 :   GEN p0, p1, q0, q1, p2, q2, P,Q, v = cgetg(3,t_VEC);
    3608             : 
    3609      609308 :   gel(v,1) = P = cgetg(lx, t_VEC);
    3610      609308 :   gel(v,2) = Q = cgetg(lx, t_VEC);
    3611      609308 :   p0 = gen_1; q0 = gen_0;
    3612      609308 :   gel(P, 1) = p1 = gel(x,1); gel(Q, 1) = q1 = gen_1;
    3613     2092209 :   for (i=2; i<lx; i++)
    3614             :   {
    3615     1482901 :     GEN a = gel(x,i);
    3616     1482901 :     gel(P, i) = p2 = addmulii(p0, a, p1); p0 = p1; p1 = p2;
    3617     1482901 :     gel(Q, i) = q2 = addmulii(q0, a, q1); q0 = q1; q1 = q2;
    3618             :   }
    3619      609308 :   return v;
    3620             : }
    3621             : 
    3622             : /* write Mod(x,N) as a/b, gcd(a,b) = 1, b <= B (no condition if B = NULL) */
    3623             : static GEN
    3624          42 : mod_to_frac(GEN x, GEN N, GEN B)
    3625             : {
    3626             :   GEN a, b, A;
    3627          42 :   if (B) A = divii(shifti(N, -1), B);
    3628             :   else
    3629             :   {
    3630          14 :     A = sqrti(shifti(N, -1));
    3631          14 :     B = A;
    3632             :   }
    3633          42 :   if (!Fp_ratlift(x, N, A,B,&a,&b) || !equali1( gcdii(a,b) )) return NULL;
    3634          28 :   return equali1(b)? a: mkfrac(a,b);
    3635             : }
    3636             : 
    3637             : static GEN
    3638          70 : mod_to_rfrac(GEN x, GEN N, long B)
    3639             : {
    3640             :   GEN a, b;
    3641          70 :   long A, d = degpol(N);
    3642          70 :   if (B >= 0) A = d-1 - B;
    3643             :   else
    3644             :   {
    3645          42 :     B = d >> 1;
    3646          42 :     A = odd(d)? B : B-1;
    3647             :   }
    3648          70 :   if (varn(N) != varn(x)) x = scalarpol(x, varn(N));
    3649          70 :   if (!RgXQ_ratlift(x, N, A, B, &a,&b) || degpol(RgX_gcd(a,b)) > 0) return NULL;
    3650          56 :   return gdiv(a,b);
    3651             : }
    3652             : 
    3653             : /* k > 0 t_INT, x a t_FRAC, returns the convergent a/b
    3654             :  * of the continued fraction of x with b <= k maximal */
    3655             : static GEN
    3656           7 : bestappr_frac(GEN x, GEN k)
    3657             : {
    3658             :   pari_sp av;
    3659             :   GEN p0, p1, p, q0, q1, q, a, y;
    3660             : 
    3661           7 :   if (cmpii(gel(x,2),k) <= 0) return gcopy(x);
    3662           0 :   av = avma; y = x;
    3663           0 :   p1 = gen_1; p0 = truedvmdii(gel(x,1), gel(x,2), &a); /* = floor(x) */
    3664           0 :   q1 = gen_0; q0 = gen_1;
    3665           0 :   x = mkfrac(a, gel(x,2)); /* = frac(x); now 0<= x < 1 */
    3666             :   for(;;)
    3667             :   {
    3668           0 :     x = ginv(x); /* > 1 */
    3669           0 :     a = typ(x)==t_INT? x: divii(gel(x,1), gel(x,2));
    3670           0 :     if (cmpii(a,k) > 0)
    3671             :     { /* next partial quotient will overflow limits */
    3672             :       GEN n, d;
    3673           0 :       a = divii(subii(k, q1), q0);
    3674           0 :       p = addii(mulii(a,p0), p1); p1=p0; p0=p;
    3675           0 :       q = addii(mulii(a,q0), q1); q1=q0; q0=q;
    3676             :       /* compare |y-p0/q0|, |y-p1/q1| */
    3677           0 :       n = gel(y,1);
    3678           0 :       d = gel(y,2);
    3679           0 :       if (abscmpii(mulii(q1, subii(mulii(q0,n), mulii(d,p0))),
    3680             :                    mulii(q0, subii(mulii(q1,n), mulii(d,p1)))) < 0)
    3681           0 :                    { p1 = p0; q1 = q0; }
    3682           0 :       break;
    3683             :     }
    3684           0 :     p = addii(mulii(a,p0), p1); p1=p0; p0=p;
    3685           0 :     q = addii(mulii(a,q0), q1); q1=q0; q0=q;
    3686             : 
    3687           0 :     if (cmpii(q0,k) > 0) break;
    3688           0 :     x = gsub(x,a); /* 0 <= x < 1 */
    3689           0 :     if (typ(x) == t_INT) { p1 = p0; q1 = q0; break; } /* x = 0 */
    3690             : 
    3691             :   }
    3692           0 :   return gerepileupto(av, gdiv(p1,q1));
    3693             : }
    3694             : /* k > 0 t_INT, x != 0 a t_REAL, returns the convergent a/b
    3695             :  * of the continued fraction of x with b <= k maximal */
    3696             : static GEN
    3697     1522524 : bestappr_real(GEN x, GEN k)
    3698             : {
    3699     1522524 :   pari_sp av = avma;
    3700     1522524 :   GEN kr, p0, p1, p, q0, q1, q, a, y = x;
    3701             : 
    3702     1522524 :   p1 = gen_1; a = p0 = floorr(x);
    3703     1522429 :   q1 = gen_0; q0 = gen_1;
    3704     1522429 :   x = subri(x,a); /* 0 <= x < 1 */
    3705     1522431 :   if (!signe(x)) { cgiv(x); return a; }
    3706     1416551 :   kr = itor(k, realprec(x));
    3707             :   for(;;)
    3708     1687472 :   {
    3709             :     long d;
    3710     3104040 :     x = invr(x); /* > 1 */
    3711     3103933 :     if (cmprr(x,kr) > 0)
    3712             :     { /* next partial quotient will overflow limits */
    3713     1403615 :       a = divii(subii(k, q1), q0);
    3714     1403551 :       p = addii(mulii(a,p0), p1); p1=p0; p0=p;
    3715     1403601 :       q = addii(mulii(a,q0), q1); q1=q0; q0=q;
    3716             :       /* compare |y-p0/q0|, |y-p1/q1| */
    3717     1403555 :       if (abscmprr(mulir(q1, subri(mulir(q0,y), p0)),
    3718             :                    mulir(q0, subri(mulir(q1,y), p1))) < 0)
    3719      124346 :                    { p1 = p0; q1 = q0; }
    3720     1403641 :       break;
    3721             :     }
    3722     1700406 :     d = nbits2prec(expo(x) + 1);
    3723     1700408 :     if (d > lg(x)) { p1 = p0; q1 = q0; break; } /* original x was ~ 0 */
    3724             : 
    3725     1700212 :     a = truncr(x); /* truncr(x) will NOT raise e_PREC */
    3726     1700174 :     p = addii(mulii(a,p0), p1); p1=p0; p0=p;
    3727     1700205 :     q = addii(mulii(a,q0), q1); q1=q0; q0=q;
    3728             : 
    3729     1700192 :     if (cmpii(q0,k) > 0) break;
    3730     1694176 :     x = subri(x,a); /* 0 <= x < 1 */
    3731     1694169 :     if (!signe(x)) { p1 = p0; q1 = q0; break; }
    3732             :   }
    3733     1416555 :   if (signe(q1) < 0) { togglesign_safe(&p1); togglesign_safe(&q1); }
    3734     1416555 :   return gerepilecopy(av, equali1(q1)? p1: mkfrac(p1,q1));
    3735             : }
    3736             : 
    3737             : /* k t_INT or NULL */
    3738             : static GEN
    3739     2510959 : bestappr_Q(GEN x, GEN k)
    3740             : {
    3741     2510959 :   long lx, tx = typ(x), i;
    3742             :   GEN a, y;
    3743             : 
    3744     2510959 :   switch(tx)
    3745             :   {
    3746         136 :     case t_INT: return icopy(x);
    3747           7 :     case t_FRAC: return k? bestappr_frac(x, k): gcopy(x);
    3748     1773762 :     case t_REAL:
    3749     1773762 :       if (!signe(x)) return gen_0;
    3750             :       /* i <= e iff nbits2lg(e+1) > lg(x) iff floorr(x) fails */
    3751     1522502 :       i = bit_prec(x); if (i <= expo(x)) return NULL;
    3752     1522522 :       return bestappr_real(x, k? k: int2n(i));
    3753             : 
    3754          28 :     case t_INTMOD: {
    3755          28 :       pari_sp av = avma;
    3756          28 :       a = mod_to_frac(gel(x,2), gel(x,1), k); if (!a) return NULL;
    3757          21 :       return gerepilecopy(av, a);
    3758             :     }
    3759          14 :     case t_PADIC: {
    3760          14 :       pari_sp av = avma;
    3761          14 :       long v = valp(x);
    3762          14 :       a = mod_to_frac(gel(x,4), gel(x,3), k); if (!a) return NULL;
    3763           7 :       if (v) a = gmul(a, powis(gel(x,2), v));
    3764           7 :       return gerepilecopy(av, a);
    3765             :     }
    3766             : 
    3767         312 :     case t_COMPLEX: {
    3768         312 :       pari_sp av = avma;
    3769         312 :       y = cgetg(3, t_COMPLEX);
    3770         312 :       gel(y,2) = bestappr(gel(x,2), k);
    3771         312 :       gel(y,1) = bestappr(gel(x,1), k);
    3772         312 :       if (gequal0(gel(y,2))) return gerepileupto(av, gel(y,1));
    3773          91 :       return y;
    3774             :     }
    3775           0 :     case t_SER:
    3776           0 :       if (ser_isexactzero(x)) return gcopy(x);
    3777             :       /* fall through */
    3778             :     case t_POLMOD: case t_POL: case t_RFRAC:
    3779             :     case t_VEC: case t_COL: case t_MAT:
    3780      736700 :       y = cgetg_copy(x, &lx);
    3781      736745 :       if (lontyp[tx] == 1) i = 1; else { y[1] = x[1]; i = 2; }
    3782     3180968 :       for (; i<lx; i++)
    3783             :       {
    3784     2444238 :         a = bestappr_Q(gel(x,i),k); if (!a) return NULL;
    3785     2444223 :         gel(y,i) = a;
    3786             :       }
    3787      736730 :       if (tx == t_POL) return normalizepol(y);
    3788      736716 :       if (tx == t_SER) return normalizeser(y);
    3789      736716 :       return y;
    3790             :   }
    3791           0 :   pari_err_TYPE("bestappr_Q",x);
    3792             :   return NULL; /* LCOV_EXCL_LINE */
    3793             : }
    3794             : 
    3795             : static GEN
    3796          56 : bestappr_ser(GEN x, long B)
    3797             : {
    3798          56 :   long dN, v = valser(x), lx = lg(x);
    3799             :   GEN t;
    3800          56 :   x = normalizepol(ser2pol_i(x, lx));
    3801          56 :   dN = lx-2;
    3802          56 :   if (v > 0)
    3803             :   {
    3804          14 :     x = RgX_shift_shallow(x, v);
    3805          14 :     dN += v;
    3806             :   }
    3807          42 :   else if (v < 0)
    3808             :   {
    3809           7 :     if (B >= 0) B = maxss(B+v, 0);
    3810             :   }
    3811          56 :   t = mod_to_rfrac(x, pol_xn(dN, varn(x)), B);
    3812          56 :   if (!t) return NULL;
    3813          42 :   if (v < 0)
    3814             :   {
    3815             :     GEN a, b;
    3816             :     long vx;
    3817           7 :     if (typ(t) == t_POL) return RgX_mulXn(t, v);
    3818             :     /* t_RFRAC */
    3819           7 :     vx = varn(x);
    3820           7 :     a = gel(t,1);
    3821           7 :     b = gel(t,2);
    3822           7 :     v -= RgX_valrem(b, &b);
    3823           7 :     if (typ(a) == t_POL && varn(a) == vx) v += RgX_valrem(a, &a);
    3824           7 :     if (v < 0) b = RgX_shift(b, -v);
    3825           0 :     else if (v > 0) {
    3826           0 :       if (typ(a) != t_POL || varn(a) != vx) a = scalarpol_shallow(a, vx);
    3827           0 :       a = RgX_shift(a, v);
    3828             :     }
    3829           7 :     t = mkrfraccopy(a, b);
    3830             :   }
    3831          42 :   return t;
    3832             : }
    3833             : static GEN bestappr_RgX(GEN x, long B);
    3834             : /* x t_POLMOD, B >= 0 or < 0 [omit condition on B].
    3835             :  * Look for coprime t_POL a,b, deg(b)<=B, such that a/b = x */
    3836             : static GEN
    3837          77 : bestappr_RgX(GEN x, long B)
    3838             : {
    3839          77 :   long i, lx, tx = typ(x);
    3840             :   GEN y, t;
    3841          77 :   switch(tx)
    3842             :   {
    3843           0 :     case t_INT: case t_REAL: case t_INTMOD: case t_FRAC:
    3844             :     case t_COMPLEX: case t_PADIC: case t_QUAD: case t_POL:
    3845           0 :       return gcopy(x);
    3846             : 
    3847          14 :     case t_RFRAC: {
    3848          14 :       pari_sp av = avma;
    3849          14 :       if (B < 0 || degpol(gel(x,2)) <= B) return gcopy(x);
    3850           7 :       x = rfrac_to_ser_i(x, 2*B+1);
    3851           7 :       t = bestappr_ser(x, B); if (!t) return NULL;
    3852           0 :       return gerepileupto(av, t);
    3853             :     }
    3854          14 :     case t_POLMOD: {
    3855          14 :       pari_sp av = avma;
    3856          14 :       t = mod_to_rfrac(gel(x,2), gel(x,1), B); if (!t) return NULL;
    3857          14 :       return gerepileupto(av, t);
    3858             :     }
    3859          49 :     case t_SER: {
    3860          49 :       pari_sp av = avma;
    3861          49 :       t = bestappr_ser(x, B); if (!t) return NULL;
    3862          42 :       return gerepileupto(av, t);
    3863             :     }
    3864             : 
    3865           0 :     case t_VEC: case t_COL: case t_MAT:
    3866           0 :       y = cgetg_copy(x, &lx);
    3867           0 :       if (lontyp[tx] == 1) i = 1; else { y[1] = x[1]; i = 2; }
    3868           0 :       for (; i<lx; i++)
    3869             :       {
    3870           0 :         t = bestappr_RgX(gel(x,i),B); if (!t) return NULL;
    3871           0 :         gel(y,i) = t;
    3872             :       }
    3873           0 :       return y;
    3874             :   }
    3875           0 :   pari_err_TYPE("bestappr_RgX",x);
    3876             :   return NULL; /* LCOV_EXCL_LINE */
    3877             : }
    3878             : 
    3879             : /* allow k = NULL: maximal accuracy */
    3880             : GEN
    3881       66711 : bestappr(GEN x, GEN k)
    3882             : {
    3883       66711 :   pari_sp av = avma;
    3884       66711 :   if (k) { /* replace by floor(k) */
    3885       66389 :     switch(typ(k))
    3886             :     {
    3887        4404 :       case t_INT:
    3888        4404 :         break;
    3889       61985 :       case t_REAL: case t_FRAC:
    3890       61985 :         k = floor_safe(k); /* left on stack for efficiency */
    3891       61986 :         if (!signe(k)) k = gen_1;
    3892       61986 :         break;
    3893           0 :       default:
    3894           0 :         pari_err_TYPE("bestappr [bound type]", k);
    3895           0 :         break;
    3896             :     }
    3897         322 :   }
    3898       66712 :   x = bestappr_Q(x, k);
    3899       66712 :   if (!x) { set_avma(av); return cgetg(1,t_VEC); }
    3900       66697 :   return x;
    3901             : }
    3902             : GEN
    3903          77 : bestapprPade(GEN x, long B)
    3904             : {
    3905          77 :   pari_sp av = avma;
    3906          77 :   GEN t = bestappr_RgX(x, B);
    3907          77 :   if (!t) { set_avma(av); return cgetg(1,t_VEC); }
    3908          63 :   return t;
    3909             : }

Generated by: LCOV version 1.14