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.14.0 lcov report (development 27775-aca467eab2) Lines: 1979 2141 92.4 %
Date: 2022-07-03 07:33:15 Functions: 206 217 94.9 %
Legend: Lines: hit not hit

          Line data    Source code
       1             : /* Copyright (C) 2000  The PARI group.
       2             : 
       3             : This file is part of the PARI/GP package.
       4             : 
       5             : PARI/GP is free software; you can redistribute it and/or modify it under the
       6             : terms of the GNU General Public License as published by the Free Software
       7             : Foundation; either version 2 of the License, or (at your option) any later
       8             : version. It is distributed in the hope that it will be useful, but WITHOUT
       9             : ANY WARRANTY WHATSOEVER.
      10             : 
      11             : Check the License for details. You should have received a copy of it, along
      12             : with the package; see the file 'COPYING'. If not, write to the Free Software
      13             : Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. */
      14             : 
      15             : /*********************************************************************/
      16             : /**                     ARITHMETIC FUNCTIONS                        **/
      17             : /**                         (first part)                            **/
      18             : /*********************************************************************/
      19             : #include "pari.h"
      20             : #include "paripriv.h"
      21             : 
      22             : #define DEBUGLEVEL DEBUGLEVEL_arith
      23             : 
      24             : /******************************************************************/
      25             : /*                 GENERATOR of (Z/mZ)*                           */
      26             : /******************************************************************/
      27             : static GEN
      28        1812 : remove2(GEN q) { long v = vali(q); return v? shifti(q, -v): q; }
      29             : static ulong
      30      132320 : 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      132319 : 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        2079 : is_gener_expo(GEN p, GEN L0)
      39             : {
      40        2079 :   GEN L, q = shifti(p,-1);
      41             :   long i, l;
      42        2079 :   if (L0) {
      43         304 :     l = lg(L0);
      44         304 :     L = cgetg(l, t_VEC);
      45             :   } else {
      46        1775 :     L0 = L = odd_prime_divisors(q);
      47        1775 :     l = lg(L);
      48             :   }
      49        9217 :   for (i=1; i<l; i++) gel(L,i) = diviiexact(q, gel(L0,i));
      50        2079 :   return L;
      51             : }
      52             : static GEN
      53      204536 : u_is_gener_expo(ulong p, GEN L0)
      54             : {
      55      204536 :   const ulong q = p >> 1;
      56             :   long i;
      57             :   GEN L;
      58      204536 :   if (!L0) L0 = u_odd_prime_divisors(q);
      59      204533 :   L = cgetg_copy(L0,&i);
      60      417954 :   while (--i) L[i] = q / uel(L0,i);
      61      204532 :   return L;
      62             : }
      63             : 
      64             : int
      65      578933 : is_gener_Fl(ulong x, ulong p, ulong p_1, GEN L)
      66             : {
      67             :   long i;
      68      578933 :   if (krouu(x, p) >= 0) return 0;
      69      525777 :   for (i=lg(L)-1; i; i--)
      70             :   {
      71      316516 :     ulong t = Fl_powu(x, uel(L,i), p);
      72      316513 :     if (t == p_1 || t == 1) return 0;
      73             :   }
      74      209261 :   return 1;
      75             : }
      76             : /* assume p prime */
      77             : ulong
      78      493895 : pgener_Fl_local(ulong p, GEN L0)
      79             : {
      80      493895 :   const pari_sp av = avma;
      81      493895 :   const ulong p_1 = p-1;
      82             :   long x;
      83             :   GEN L;
      84      493895 :   if (p <= 19) switch(p)
      85             :   { /* quick trivial cases */
      86          21 :     case 2:  return 1;
      87       60564 :     case 7:
      88       60564 :     case 17: return 3;
      89      228837 :     default: return 2;
      90             :   }
      91      204473 :   L = u_is_gener_expo(p,L0);
      92      204501 :   for (x = 2;; x++)
      93      571093 :     if (is_gener_Fl(x,p,p_1,L)) return gc_ulong(av, x);
      94             : }
      95             : ulong
      96      201834 : 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        8472 : is_gener_Fp(GEN x, GEN p, GEN p_1, GEN L)
     102             : {
     103        8472 :   long i, t = lgefint(x)==3? kroui(x[2], p): kronecker(x, p);
     104        8472 :   if (t >= 0) return 0;
     105       16582 :   for (i = lg(L)-1; i; i--)
     106             :   {
     107       11904 :     GEN t = Fp_pow(x, gel(L,i), p);
     108       11904 :     if (equalii(t, p_1) || equali1(t)) return 0;
     109             :   }
     110        4678 :   return 1;
     111             : }
     112             : 
     113             : /* assume p prime, return a generator of all L[i]-Sylows in F_p^*. */
     114             : GEN
     115      354672 : pgener_Fp_local(GEN p, GEN L0)
     116             : {
     117      354672 :   pari_sp av0 = avma;
     118             :   GEN x, p_1, L;
     119      354672 :   if (lgefint(p) == 3)
     120             :   {
     121             :     ulong z;
     122      352600 :     if (p[2] == 2) return gen_1;
     123      257778 :     if (L0) L0 = ZV_to_nv(L0);
     124      257771 :     z = pgener_Fl_local(uel(p,2), L0);
     125      257819 :     set_avma(av0); return utoipos(z);
     126             :   }
     127        2072 :   p_1 = subiu(p,1); L = is_gener_expo(p, L0);
     128        2074 :   x = utoipos(2);
     129        4443 :   for (;; x[2]++) { if (is_gener_Fp(x, p, p_1, L)) break; }
     130        2074 :   set_avma(av0); return utoipos(uel(x,2));
     131             : }
     132             : 
     133             : GEN
     134       44016 : pgener_Fp(GEN p) { return pgener_Fp_local(p, NULL); }
     135             : 
     136             : ulong
     137      114214 : pgener_Zl(ulong p)
     138             : {
     139      114214 :   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      114214 :   if (p == 40487) return 10;
     142             : #ifndef LONG_IS_64BIT
     143       16312 :   return pgener_Fl(p);
     144             : #else
     145       97902 :   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          30 :     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      113687 : pgener_Zp(GEN p)
     162             : {
     163      113687 :   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 :     set_avma(av); return utoipos(uel(x,2));
     172             :   }
     173             : }
     174             : 
     175             : static GEN
     176         245 : gener_Zp(GEN q, GEN F)
     177             : {
     178         245 :   GEN p = NULL;
     179         245 :   long e = 0;
     180         245 :   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", "argument","=",F,F);
     189          14 :       e = itos(gel(E,i));
     190             :     }
     191          14 :     if (!p) pari_err_DOMAIN("znprimroot", "argument","=",F,F);
     192             :   }
     193             :   else
     194         231 :     e = Z_isanypower(q, &p);
     195         245 :   return e > 1? pgener_Zp(p): pgener_Fp(q);
     196             : }
     197             : 
     198             : GEN
     199         315 : znprimroot(GEN N)
     200             : {
     201         315 :   pari_sp av = avma;
     202             :   GEN x, n, F;
     203             : 
     204         315 :   if ((F = check_arith_non0(N,"znprimroot")))
     205             :   {
     206          14 :     F = clean_Z_factor(F);
     207          14 :     N = typ(N) == t_VEC? gel(N,1): factorback(F);
     208             :   }
     209         308 :   N = absi_shallow(N);
     210         308 :   if (abscmpiu(N, 4) <= 0) { set_avma(av); return mkintmodu(N[2]-1,N[2]); }
     211         259 :   switch(mod4(N))
     212             :   {
     213          14 :     case 0: /* N = 0 mod 4 */
     214          14 :       pari_err_DOMAIN("znprimroot", "argument","=",N,N);
     215           0 :       x = NULL; break;
     216          21 :     case 2: /* N = 2 mod 4 */
     217          21 :       n = shifti(N,-1); /* becomes odd */
     218          21 :       x = gener_Zp(n,F); if (!mod2(x)) x = addii(x,n);
     219          21 :       break;
     220         224 :     default: /* N odd */
     221         224 :       x = gener_Zp(N,F);
     222         224 :       break;
     223             :   }
     224         245 :   return gerepilecopy(av, mkintmod(x, N));
     225             : }
     226             : 
     227             : /* n | (p-1), returns a primitive n-th root of 1 in F_p^* */
     228             : GEN
     229           0 : rootsof1_Fp(GEN n, GEN p)
     230             : {
     231           0 :   pari_sp av = avma;
     232           0 :   GEN L = odd_prime_divisors(n); /* 2 implicit in pgener_Fp_local */
     233           0 :   GEN z = pgener_Fp_local(p, L);
     234           0 :   z = Fp_pow(z, diviiexact(subiu(p,1), n), p); /* prim. n-th root of 1 */
     235           0 :   return gerepileuptoint(av, z);
     236             : }
     237             : 
     238             : GEN
     239         217 : rootsof1u_Fp(ulong n, GEN p)
     240             : {
     241         217 :   pari_sp av = avma;
     242         217 :   GEN z, L = u_odd_prime_divisors(n); /* 2 implicit in pgener_Fp_local */
     243         217 :   z = pgener_Fp_local(p, Flv_to_ZV(L));
     244         217 :   z = Fp_pow(z, diviuexact(subiu(p,1), n), p); /* prim. n-th root of 1 */
     245         217 :   return gerepileuptoint(av, z);
     246             : }
     247             : 
     248             : ulong
     249       32115 : rootsof1_Fl(ulong n, ulong p)
     250             : {
     251       32115 :   pari_sp av = avma;
     252       32115 :   GEN L = u_odd_prime_divisors(n); /* 2 implicit in pgener_Fl_local */
     253       32115 :   ulong z = pgener_Fl_local(p, L);
     254       32115 :   z = Fl_powu(z, (p-1) / n, p); /* prim. n-th root of 1 */
     255       32115 :   return gc_ulong(av,z);
     256             : }
     257             : 
     258             : /*********************************************************************/
     259             : /**                     INVERSE TOTIENT FUNCTION                    **/
     260             : /*********************************************************************/
     261             : /* N t_INT, L a ZV containing all prime divisors of N, and possibly other
     262             :  * primes. Return factor(N) */
     263             : GEN
     264      350651 : Z_factor_listP(GEN N, GEN L)
     265             : {
     266      350651 :   long i, k, l = lg(L);
     267      350651 :   GEN P = cgetg(l, t_COL), E = cgetg(l, t_COL);
     268     1346688 :   for (i = k = 1; i < l; i++)
     269             :   {
     270      996037 :     GEN p = gel(L,i);
     271      996037 :     long v = Z_pvalrem(N, p, &N);
     272      996037 :     if (v)
     273             :     {
     274      792176 :       gel(P,k) = p;
     275      792176 :       gel(E,k) = utoipos(v);
     276      792176 :       k++;
     277             :     }
     278             :   }
     279      350651 :   setlg(P, k);
     280      350651 :   setlg(E, k); return mkmat2(P,E);
     281             : }
     282             : 
     283             : /* look for x such that phi(x) = n, p | x => p > m (if m = NULL: no condition).
     284             :  * L is a list of primes containing all prime divisors of n. */
     285             : static long
     286      621565 : istotient_i(GEN n, GEN m, GEN L, GEN *px)
     287             : {
     288      621565 :   pari_sp av = avma, av2;
     289             :   GEN k, D;
     290             :   long i, v;
     291      621565 :   if (m && mod2(n))
     292             :   {
     293      270914 :     if (!equali1(n)) return 0;
     294       69986 :     if (px) *px = gen_1;
     295       69986 :     return 1;
     296             :   }
     297      350651 :   D = divisors(Z_factor_listP(shifti(n, -1), L));
     298             :   /* loop through primes p > m, d = p-1 | n */
     299      350651 :   av2 = avma;
     300      350651 :   if (!m)
     301             :   { /* special case p = 2, d = 1 */
     302       69986 :     k = n;
     303       69986 :     for (v = 1;; v++) {
     304       69986 :       if (istotient_i(k, gen_2, L, px)) {
     305       69986 :         if (px) *px = shifti(*px, v);
     306       69986 :         return 1;
     307             :       }
     308           0 :       if (mod2(k)) break;
     309           0 :       k = shifti(k,-1);
     310             :     }
     311           0 :     set_avma(av2);
     312             :   }
     313     1099462 :   for (i = 1; i < lg(D); ++i)
     314             :   {
     315     1001588 :     GEN p, d = shifti(gel(D, i), 1); /* even divisors of n */
     316     1001588 :     if (m && cmpii(d, m) < 0) continue;
     317      677782 :     p = addiu(d, 1);
     318      677782 :     if (!isprime(p)) continue;
     319      442064 :     k = diviiexact(n, d);
     320      481593 :     for (v = 1;; v++) {
     321             :       GEN r;
     322      481593 :       if (istotient_i(k, p, L, px)) {
     323      182791 :         if (px) *px = mulii(*px, powiu(p, v));
     324      182791 :         return 1;
     325             :       }
     326      298802 :       k = dvmdii(k, p, &r);
     327      298802 :       if (r != gen_0) break;
     328             :     }
     329      259273 :     set_avma(av2);
     330             :   }
     331       97874 :   return gc_long(av,0);
     332             : }
     333             : 
     334             : /* find x such that phi(x) = n */
     335             : long
     336       70000 : istotient(GEN n, GEN *px)
     337             : {
     338       70000 :   pari_sp av = avma;
     339       70000 :   if (typ(n) != t_INT) pari_err_TYPE("istotient", n);
     340       70000 :   if (signe(n) < 1) return 0;
     341       70000 :   if (mod2(n))
     342             :   {
     343          14 :     if (!equali1(n)) return 0;
     344          14 :     if (px) *px = gen_1;
     345          14 :     return 1;
     346             :   }
     347       69986 :   if (istotient_i(n, NULL, gel(Z_factor(n), 1), px))
     348             :   {
     349       69986 :     if (!px) set_avma(av);
     350             :     else
     351       69986 :       *px = gerepileuptoint(av, *px);
     352       69986 :     return 1;
     353             :   }
     354           0 :   return gc_long(av,0);
     355             : }
     356             : 
     357             : /*********************************************************************/
     358             : /**                        KRONECKER SYMBOL                         **/
     359             : /*********************************************************************/
     360             : /* t = 3,5 mod 8 ?  (= 2 not a square mod t) */
     361             : static int
     362   298297947 : ome(long t)
     363             : {
     364   298297947 :   switch(t & 7)
     365             :   {
     366   169433605 :     case 3:
     367   169433605 :     case 5: return 1;
     368   128864342 :     default: return 0;
     369             :   }
     370             : }
     371             : /* t a t_INT, is t = 3,5 mod 8 ? */
     372             : static int
     373     4098386 : gome(GEN t)
     374     4098386 : { return signe(t)? ome( mod2BIL(t) ): 0; }
     375             : 
     376             : /* assume y odd, return kronecker(x,y) * s */
     377             : static long
     378   212697992 : krouu_s(ulong x, ulong y, long s)
     379             : {
     380   212697992 :   ulong x1 = x, y1 = y, z;
     381   955723329 :   while (x1)
     382             :   {
     383   743045056 :     long r = vals(x1);
     384   743006678 :     if (r)
     385             :     {
     386   395412975 :       if (odd(r) && ome(y1)) s = -s;
     387   395431634 :       x1 >>= r;
     388             :     }
     389   743025337 :     if (x1 & y1 & 2) s = -s;
     390   743025337 :     z = y1 % x1; y1 = x1; x1 = z;
     391             :   }
     392   212678273 :   return (y1 == 1)? s: 0;
     393             : }
     394             : 
     395             : long
     396     9518524 : kronecker(GEN x, GEN y)
     397             : {
     398     9518524 :   pari_sp av = avma;
     399     9518524 :   long s = 1, r;
     400             :   ulong xu;
     401             : 
     402     9518524 :   if (typ(x) != t_INT) pari_err_TYPE("kronecker",x);
     403     9518524 :   if (typ(y) != t_INT) pari_err_TYPE("kronecker",y);
     404     9518524 :   switch (signe(y))
     405             :   {
     406          63 :     case -1: y = negi(y); if (signe(x) < 0) s = -1; break;
     407           0 :     case 0: return is_pm1(x);
     408             :   }
     409     9518524 :   r = vali(y);
     410     9518526 :   if (r)
     411             :   {
     412      312749 :     if (!mpodd(x)) return gc_long(av,0);
     413      310887 :     if (odd(r) && gome(x)) s = -s;
     414      310887 :     y = shifti(y,-r);
     415             :   }
     416     9516663 :   x = modii(x,y);
     417    10517635 :   while (lgefint(x) > 3) /* x < y */
     418             :   {
     419             :     GEN z;
     420     1001025 :     r = vali(x);
     421     1001002 :     if (r)
     422             :     {
     423      549381 :       if (odd(r) && gome(y)) s = -s;
     424      549345 :       x = shifti(x,-r);
     425             :     }
     426             :     /* x=3 mod 4 && y=3 mod 4 ? (both are odd here) */
     427      999655 :     if (mod2BIL(x) & mod2BIL(y) & 2) s = -s;
     428     1000008 :     z = remii(y,x); y = x; x = z;
     429     1001124 :     if (gc_needed(av,2))
     430             :     {
     431           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"kronecker");
     432           0 :       gerepileall(av, 2, &x, &y);
     433             :     }
     434             :   }
     435     9516610 :   xu = itou(x);
     436     9516613 :   if (!xu) return is_pm1(y)? s: 0;
     437     9465082 :   r = vals(xu);
     438     9465095 :   if (r)
     439             :   {
     440     4705412 :     if (odd(r) && gome(y)) s = -s;
     441     4705412 :     xu >>= r;
     442             :   }
     443             :   /* x=3 mod 4 && y=3 mod 4 ? (both are odd here) */
     444     9465095 :   if (xu & mod2BIL(y) & 2) s = -s;
     445     9465107 :   return gc_long(av, krouu_s(umodiu(y,xu), xu, s));
     446             : }
     447             : 
     448             : long
     449       35973 : krois(GEN x, long y)
     450             : {
     451             :   ulong yu;
     452       35973 :   long s = 1;
     453             : 
     454       35973 :   if (y <= 0)
     455             :   {
     456           7 :     if (y == 0) return is_pm1(x);
     457           0 :     yu = (ulong)-y; if (signe(x) < 0) s = -1;
     458             :   }
     459             :   else
     460       35966 :     yu = (ulong)y;
     461       35966 :   if (!odd(yu))
     462             :   {
     463             :     long r;
     464       14938 :     if (!mpodd(x)) return 0;
     465       11088 :     r = vals(yu); yu >>= r;
     466       11088 :     if (odd(r) && gome(x)) s = -s;
     467             :   }
     468       32116 :   return krouu_s(umodiu(x, yu), yu, s);
     469             : }
     470             : /* assume y != 0 */
     471             : long
     472    27339887 : kroiu(GEN x, ulong y)
     473             : {
     474             :   long r;
     475    27339887 :   if (odd(y)) return krouu_s(umodiu(x,y), y, 1);
     476      297091 :   if (!mpodd(x)) return 0;
     477      203131 :   r = vals(y); y >>= r;
     478      203134 :   return krouu_s(umodiu(x,y), y, (odd(r) && gome(x))? -1: 1);
     479             : }
     480             : 
     481             : /* assume y > 0, odd, return s * kronecker(x,y) */
     482             : static long
     483      163172 : krouodd(ulong x, GEN y, long s)
     484             : {
     485             :   long r;
     486      163172 :   if (lgefint(y) == 3) return krouu_s(x, y[2], s);
     487       19820 :   if (!x) return 0; /* y != 1 */
     488       19820 :   r = vals(x);
     489       19820 :   if (r)
     490             :   {
     491       10220 :     if (odd(r) && gome(y)) s = -s;
     492       10220 :     x >>= r;
     493             :   }
     494             :   /* x=3 mod 4 && y=3 mod 4 ? (both are odd here) */
     495       19820 :   if (x & mod2BIL(y) & 2) s = -s;
     496       19820 :   return krouu_s(umodiu(y,x), x, s);
     497             : }
     498             : 
     499             : long
     500      154767 : krosi(long x, GEN y)
     501             : {
     502      154767 :   const pari_sp av = avma;
     503      154767 :   long s = 1, r;
     504      154767 :   switch (signe(y))
     505             :   {
     506           0 :     case -1: y = negi(y); if (x < 0) s = -1; break;
     507           0 :     case 0: return (x==1 || x==-1);
     508             :   }
     509      154767 :   r = vali(y);
     510      154767 :   if (r)
     511             :   {
     512       16884 :     if (!odd(x)) return gc_long(av,0);
     513       16884 :     if (odd(r) && ome(x)) s = -s;
     514       16884 :     y = shifti(y,-r);
     515             :   }
     516      154767 :   if (x < 0) { x = -x; if (mod4(y) == 3) s = -s; }
     517      154767 :   return gc_long(av, krouodd((ulong)x, y, s));
     518             : }
     519             : 
     520             : long
     521        8405 : kroui(ulong x, GEN y)
     522             : {
     523        8405 :   const pari_sp av = avma;
     524        8405 :   long s = 1, r;
     525        8405 :   switch (signe(y))
     526             :   {
     527           0 :     case -1: y = negi(y); break;
     528           0 :     case 0: return x==1UL;
     529             :   }
     530        8405 :   r = vali(y);
     531        8405 :   if (r)
     532             :   {
     533           0 :     if (!odd(x)) return gc_long(av,0);
     534           0 :     if (odd(r) && ome(x)) s = -s;
     535           0 :     y = shifti(y,-r);
     536             :   }
     537        8405 :   return gc_long(av,  krouodd(x, y, s));
     538             : }
     539             : 
     540             : long
     541    94633498 : kross(long x, long y)
     542             : {
     543             :   ulong yu;
     544    94633498 :   long s = 1;
     545             : 
     546    94633498 :   if (y <= 0)
     547             :   {
     548       68943 :     if (y == 0) return (labs(x)==1);
     549       68915 :     yu = (ulong)-y; if (x < 0) s = -1;
     550             :   }
     551             :   else
     552    94564555 :     yu = (ulong)y;
     553    94633470 :   if (!odd(yu))
     554             :   {
     555             :     long r;
     556    24360810 :     if (!odd(x)) return 0;
     557    17450540 :     r = vals(yu); yu >>= r;
     558    17450540 :     if (odd(r) && ome(x)) s = -s;
     559             :   }
     560    87723200 :   x %= (long)yu; if (x < 0) x += yu;
     561    87723200 :   return krouu_s((ulong)x, yu, s);
     562             : }
     563             : 
     564             : long
     565    88084008 : krouu(ulong x, ulong y)
     566             : {
     567             :   long r;
     568    88084008 :   if (odd(y)) return krouu_s(x, y, 1);
     569        5707 :   if (!odd(x)) return 0;
     570        6182 :   r = vals(y); y >>= r;
     571        6182 :   return krouu_s(x, y, (odd(r) && ome(x))? -1: 1);
     572             : }
     573             : 
     574             : /*********************************************************************/
     575             : /**                          HILBERT SYMBOL                         **/
     576             : /*********************************************************************/
     577             : /* x,y are t_INT or t_REAL */
     578             : static long
     579        7343 : mphilbertoo(GEN x, GEN y)
     580             : {
     581        7343 :   long sx = signe(x), sy = signe(y);
     582        7343 :   if (!sx || !sy) return 0;
     583        7343 :   return (sx < 0 && sy < 0)? -1: 1;
     584             : }
     585             : 
     586             : long
     587      134624 : hilbertii(GEN x, GEN y, GEN p)
     588             : {
     589             :   pari_sp av;
     590             :   long oddvx, oddvy, z;
     591             : 
     592      134624 :   if (!p) return mphilbertoo(x,y);
     593      127302 :   if (is_pm1(p) || signe(p) < 0) pari_err_PRIME("hilbertii",p);
     594      127302 :   if (!signe(x) || !signe(y)) return 0;
     595      127281 :   av = avma;
     596      127281 :   oddvx = odd(Z_pvalrem(x,p,&x));
     597      127281 :   oddvy = odd(Z_pvalrem(y,p,&y));
     598             :   /* x, y are p-units, compute hilbert(x * p^oddvx, y * p^oddvy, p) */
     599      127281 :   if (absequaliu(p, 2))
     600             :   {
     601       12019 :     z = (Mod4(x) == 3 && Mod4(y) == 3)? -1: 1;
     602       12019 :     if (oddvx && gome(y)) z = -z;
     603       12019 :     if (oddvy && gome(x)) z = -z;
     604             :   }
     605             :   else
     606             :   {
     607      115262 :     z = (oddvx && oddvy && mod4(p) == 3)? -1: 1;
     608      115262 :     if (oddvx && kronecker(y,p) < 0) z = -z;
     609      115262 :     if (oddvy && kronecker(x,p) < 0) z = -z;
     610             :   }
     611      127281 :   return gc_long(av, z);
     612             : }
     613             : 
     614             : static void
     615         196 : err_prec(void) { pari_err_PREC("hilbert"); }
     616             : static void
     617         161 : err_p(GEN p, GEN q) { pari_err_MODULUS("hilbert", p,q); }
     618             : static void
     619          56 : err_oo(GEN p) { pari_err_MODULUS("hilbert", p, strtoGENstr("oo")); }
     620             : 
     621             : /* x t_INTMOD, *pp = prime or NULL [ unset, set it to x.mod ].
     622             :  * Return lift(x) provided it's p-adic accuracy is large enough to decide
     623             :  * hilbert()'s value [ problem at p = 2 ] */
     624             : static GEN
     625         420 : lift_intmod(GEN x, GEN *pp)
     626             : {
     627         420 :   GEN p = *pp, N = gel(x,1);
     628         420 :   x = gel(x,2);
     629         420 :   if (!p)
     630             :   {
     631         266 :     *pp = p = N;
     632         266 :     switch(itos_or_0(p))
     633             :     {
     634         126 :       case 2:
     635         126 :       case 4: err_prec();
     636             :     }
     637         140 :     return x;
     638             :   }
     639         154 :   if (!signe(p)) err_oo(N);
     640         112 :   if (absequaliu(p,2))
     641          42 :   { if (vali(N) <= 2) err_prec(); }
     642             :   else
     643          70 :   { if (!dvdii(N,p)) err_p(N,p); }
     644          28 :   if (!signe(x)) err_prec();
     645          21 :   return x;
     646             : }
     647             : /* x t_PADIC, *pp = prime or NULL [ unset, set it to x.p ].
     648             :  * Return lift(x)*p^(v(x) mod 2) provided it's p-adic accuracy is large enough
     649             :  * to decide hilbert()'s value [ problem at p = 2 ]*/
     650             : static GEN
     651         210 : lift_padic(GEN x, GEN *pp)
     652             : {
     653         210 :   GEN p = *pp, q = gel(x,2), y = gel(x,4);
     654         210 :   if (!p) *pp = p = q;
     655         147 :   else if (!equalii(p,q)) err_p(p, q);
     656         105 :   if (absequaliu(p,2) && precp(x) <= 2) err_prec();
     657          70 :   if (!signe(y)) err_prec();
     658          70 :   return odd(valp(x))? mulii(p,y): y;
     659             : }
     660             : 
     661             : long
     662       59857 : hilbert(GEN x, GEN y, GEN p)
     663             : {
     664       59857 :   pari_sp av = avma;
     665       59857 :   long tx = typ(x), ty = typ(y);
     666             : 
     667       59857 :   if (p && typ(p) != t_INT) pari_err_TYPE("hilbert",p);
     668       59857 :   if (tx == t_REAL)
     669             :   {
     670          77 :     if (p && signe(p)) err_oo(p);
     671          63 :     switch (ty)
     672             :     {
     673           7 :       case t_INT:
     674           7 :       case t_REAL: return mphilbertoo(x,y);
     675           0 :       case t_FRAC: return mphilbertoo(x,gel(y,1));
     676          56 :       default: pari_err_TYPE2("hilbert",x,y);
     677             :     }
     678             :   }
     679       59780 :   if (ty == t_REAL)
     680             :   {
     681          14 :     if (p && signe(p)) err_oo(p);
     682          14 :     switch (tx)
     683             :     {
     684          14 :       case t_INT:
     685          14 :       case t_REAL: return mphilbertoo(x,y);
     686           0 :       case t_FRAC: return mphilbertoo(gel(x,1),y);
     687           0 :       default: pari_err_TYPE2("hilbert",x,y);
     688             :     }
     689             :   }
     690       59766 :   if (tx == t_INTMOD) { x = lift_intmod(x, &p); tx = t_INT; }
     691       59563 :   if (ty == t_INTMOD) { y = lift_intmod(y, &p); ty = t_INT; }
     692             : 
     693       59507 :   if (tx == t_PADIC) { x = lift_padic(x, &p); tx = t_INT; }
     694       59444 :   if (ty == t_PADIC) { y = lift_padic(y, &p); ty = t_INT; }
     695             : 
     696       59367 :   if (tx == t_FRAC) { tx = t_INT; x = p? mulii(gel(x,1),gel(x,2)): gel(x,1); }
     697       59367 :   if (ty == t_FRAC) { ty = t_INT; y = p? mulii(gel(y,1),gel(y,2)): gel(y,1); }
     698             : 
     699       59367 :   if (tx != t_INT || ty != t_INT) pari_err_TYPE2("hilbert",x,y);
     700       59367 :   if (p && !signe(p)) p = NULL;
     701       59367 :   return gc_long(av, hilbertii(x,y,p));
     702             : }
     703             : 
     704             : /*******************************************************************/
     705             : /*                       SQUARE ROOT MODULO p                      */
     706             : /*******************************************************************/
     707             : static void
     708     1743913 : checkp(ulong q, ulong p)
     709     1743913 : { if (!q) pari_err_PRIME("Fl_nonsquare",utoipos(p)); }
     710             : /* p = 1 (mod 4) prime, return the first quadratic nonresidue, a prime */
     711             : static ulong
     712    10368671 : nonsquare1_Fl(ulong p)
     713             : {
     714             :   forprime_t S;
     715             :   ulong q;
     716    10368671 :   if ((p & 7UL) != 1) return 2UL;
     717     3745626 :   q = p % 3; if (q == 2) return 3UL;
     718     1127340 :   checkp(q, p);
     719     1127412 :   q = p % 5; if (q == 2 || q == 3) return 5UL;
     720      396420 :   checkp(q, p);
     721      396420 :   q = p % 7; if (q != 4 && q >= 3) return 7UL;
     722      138826 :   checkp(q, p);
     723             :   /* log^2(2^64) < 1968 is enough under GRH (and p^(1/4)log(p) without it)*/
     724      138826 :   u_forprime_init(&S, 11, 1967);
     725      220074 :   while ((q = u_forprime_next(&S)))
     726             :   {
     727      220074 :     if (krouu(q, p) < 0) return q;
     728       81248 :     checkp(q, p);
     729             :   }
     730           0 :   checkp(0, p);
     731             :   return 0; /*LCOV_EXCL_LINE*/
     732             : }
     733             : /* p > 2 a prime */
     734             : ulong
     735        7853 : nonsquare_Fl(ulong p)
     736        7853 : { return ((p & 3UL) == 3)? p-1: nonsquare1_Fl(p); }
     737             : 
     738             : ulong
     739      155880 : Fl_2gener_pre(ulong p, ulong pi)
     740             : {
     741      155880 :   ulong p1 = p-1;
     742      155880 :   long e = vals(p1);
     743      155879 :   if (e == 1) return p1;
     744       58337 :   return Fl_powu_pre(nonsquare1_Fl(p), p1 >> e, p, pi);
     745             : }
     746             : 
     747             : /* Tonelli-Shanks. Assume p is prime and (a,p) != -1. */
     748             : ulong
     749    30272523 : Fl_sqrt_pre_i(ulong a, ulong y, ulong p, ulong pi)
     750             : {
     751             :   long i, e, k;
     752             :   ulong p1, q, v, w;
     753             : 
     754    30272523 :   if (!a) return 0;
     755    28930522 :   p1 = p - 1; e = vals(p1);
     756    28939656 :   if (e == 0) /* p = 2 */
     757             :   {
     758      420476 :     if (p != 2) pari_err_PRIME("Fl_sqrt [modulus]",utoi(p));
     759      421973 :     return ((a & 1) == 0)? 0: 1;
     760             :   }
     761    28519180 :   if (e == 1)
     762             :   {
     763    18207514 :     v = Fl_powu_pre(a, (p+1) >> 2, p, pi);
     764    18168092 :     if (Fl_sqr_pre(v, p, pi) != a) return ~0UL;
     765    18175483 :     p1 = p - v; if (v > p1) v = p1;
     766    18175483 :     return v;
     767             :   }
     768    10311666 :   q = p1 >> e; /* q = (p-1)/2^oo is odd */
     769    10311666 :   p1 = Fl_powu_pre(a, q >> 1, p, pi); /* a ^ [(q-1)/2] */
     770    10312036 :   if (!p1) return 0;
     771    10312036 :   v = Fl_mul_pre(a, p1, p, pi);
     772    10312124 :   w = Fl_mul_pre(v, p1, p, pi);
     773    10312140 :   if (!y) y = Fl_powu_pre(nonsquare1_Fl(p), q, p, pi);
     774    18574622 :   while (w != 1)
     775             :   { /* a*w = v^2, y primitive 2^e-th root of 1
     776             :        a square --> w even power of y, hence w^(2^(e-1)) = 1 */
     777     8264460 :     p1 = Fl_sqr_pre(w,p,pi);
     778    14746911 :     for (k=1; p1 != 1 && k < e; k++) p1 = Fl_sqr_pre(p1,p,pi);
     779     8264507 :     if (k == e) return ~0UL;
     780             :     /* w ^ (2^k) = 1 --> w = y ^ (u * 2^(e-k)), u odd */
     781     8262547 :     p1 = y;
     782    10856836 :     for (i=1; i < e-k; i++) p1 = Fl_sqr_pre(p1, p, pi);
     783     8262546 :     y = Fl_sqr_pre(p1, p, pi); e = k;
     784     8262556 :     w = Fl_mul_pre(y, w, p, pi);
     785     8262554 :     v = Fl_mul_pre(v, p1, p, pi);
     786             :   }
     787    10310162 :   p1 = p - v; if (v > p1) v = p1;
     788    10310162 :   return v;
     789             : }
     790             : 
     791             : ulong
     792    26769307 : Fl_sqrt(ulong a, ulong p)
     793    26769307 : { ulong pi = get_Fl_red(p); return Fl_sqrt_pre_i(a, 0, p, pi); }
     794             : 
     795             : ulong
     796     3453453 : Fl_sqrt_pre(ulong a, ulong p, ulong pi)
     797     3453453 : { return Fl_sqrt_pre_i(a, 0, p, pi); }
     798             : 
     799             : static ulong
     800       48606 : Fl_lgener_pre_all(ulong l, long e, ulong r, ulong p, ulong pi, ulong *pt_m)
     801             : {
     802             :   ulong x, y, m;
     803       48606 :   ulong le1 = upowuu(l, e-1);
     804       48606 :   for (x = 2; ; x++)
     805             :   {
     806       76507 :     y = Fl_powu_pre(x, r, p, pi);
     807       76508 :     if (y==1) continue;
     808       59487 :     m = Fl_powu_pre(y, le1, p, pi);
     809       59486 :     if (m != 1) break;
     810             :   }
     811       48606 :   *pt_m = m; return y;
     812             : }
     813             : 
     814             : /* solve x^l = a , l prime in G of order q.
     815             :  *
     816             :  * q =  (l^e)*r, e >= 1, (r,l) = 1
     817             :  * y generates the l-Sylow of G
     818             :  * m = y^(l^(e-1)) != 1 */
     819             : static ulong
     820      117063 : Fl_sqrtl_raw(ulong a, ulong l, ulong e, ulong r, ulong p, ulong pi, ulong y, ulong m)
     821             : {
     822             :   ulong p1, v, w, z, dl;
     823             :   ulong u2;
     824      117063 :   if (a==0) return a;
     825      117063 :   u2 = Fl_inv(l%r, r);
     826      117064 :   v = Fl_powu_pre(a, u2, p, pi);
     827      117065 :   w = Fl_powu_pre(v, l, p, pi);
     828      117064 :   w = Fl_mul_pre(w, Fl_inv(a, p), p, pi);
     829      117047 :   if (w==1) return v;
     830       47241 :   if (y==0) y = Fl_lgener_pre_all(l, e, r, p, pi, &m);
     831       67501 :   while (w!=1)
     832             :   {
     833       51795 :     ulong k = 0;
     834       51795 :     p1 = w;
     835             :     do
     836             :     {
     837       76829 :       z = p1; p1 = Fl_powu_pre(p1, l, p, pi);
     838       76829 :       if (++k == e) return ULONG_MAX;
     839       45294 :     } while (p1!=1);
     840       20260 :     dl = Fl_log_pre(z, m, l, p, pi);
     841       20260 :     dl = Fl_neg(dl, l);
     842       20260 :     p1 = Fl_powu_pre(y,dl*upowuu(l,e-k-1),p,pi);
     843       20260 :     m = Fl_powu_pre(m, dl, p, pi);
     844       20260 :     e = k;
     845       20260 :     v = Fl_mul_pre(p1,v,p,pi);
     846       20260 :     y = Fl_powu_pre(p1,l,p,pi);
     847       20260 :     w = Fl_mul_pre(y,w,p,pi);
     848             :   }
     849       15706 :   return v;
     850             : }
     851             : 
     852             : static ulong
     853      114713 : Fl_sqrtl_i(ulong a, ulong l, ulong p, ulong pi, ulong y, ulong m)
     854             : {
     855      114713 :   ulong r, e = u_lvalrem(p-1, l, &r);
     856      114712 :   return Fl_sqrtl_raw(a, l, e, r, p, pi, y, m);
     857             : }
     858             : 
     859             : ulong
     860      114713 : Fl_sqrtl_pre(ulong a, ulong l, ulong p, ulong pi)
     861      114713 : { return Fl_sqrtl_i(a, l, p, pi, 0, 0); }
     862             : 
     863             : ulong
     864           0 : Fl_sqrtl(ulong a, ulong l, ulong p)
     865           0 : { ulong pi = get_Fl_red(p); return Fl_sqrtl_i(a, l, p, pi, 0, 0); }
     866             : 
     867             : ulong
     868       85834 : Fl_sqrtn_pre(ulong a, long n, ulong p, ulong pi, ulong *zetan)
     869             : {
     870       85834 :   ulong m, q = p-1, z;
     871       85834 :   ulong nn = n >= 0 ? (ulong)n: -(ulong)n;
     872       85834 :   if (a==0)
     873             :   {
     874       48139 :     if (n < 0) pari_err_INV("Fl_sqrtn", mkintmod(gen_0,utoi(p)));
     875       48132 :     if (zetan) *zetan = 1UL;
     876       48132 :     return 0;
     877             :   }
     878       37695 :   if (n==1)
     879             :   {
     880         420 :     if (zetan) *zetan = 1;
     881         420 :     return n < 0? Fl_inv(a,p): a;
     882             :   }
     883       37275 :   if (n==2)
     884             :   {
     885        5887 :     if (zetan) *zetan = p-1;
     886        5887 :     return Fl_sqrt_pre_i(a, 0, p, pi);
     887             :   }
     888       31388 :   if (a == 1 && !zetan) return a;
     889        7969 :   m = ugcd(nn, q);
     890        7969 :   z = 1;
     891        7969 :   if (m!=1)
     892             :   {
     893        1396 :     GEN F = factoru(m);
     894             :     long i, j, e;
     895             :     ulong r, zeta, y, l;
     896        3118 :     for (i = nbrows(F); i; i--)
     897             :     {
     898        1785 :       l = ucoeff(F,i,1);
     899        1785 :       j = ucoeff(F,i,2);
     900        1785 :       e = u_lvalrem(q,l, &r);
     901        1785 :       y = Fl_lgener_pre_all(l, e, r, p, pi, &zeta);
     902        1785 :       if (zetan)
     903          98 :         z = Fl_mul_pre(z, Fl_powu_pre(y, upowuu(l,e-j), p, pi), p, pi);
     904        1785 :       if (a!=1)
     905             :         do
     906             :         {
     907        2351 :           a = Fl_sqrtl_raw(a, l, e, r, p, pi, y, zeta);
     908        2337 :           if (a==ULONG_MAX) return ULONG_MAX;
     909        2288 :         } while (--j);
     910             :     }
     911             :   }
     912        7906 :   if (m != nn)
     913             :   {
     914        6594 :     ulong qm = q/m, nm = nn/m;
     915        6594 :     a = Fl_powu_pre(a, Fl_inv(nm%qm, qm), p, pi);
     916             :   }
     917        7906 :   if (n < 0) a = Fl_inv(a, p);
     918        7906 :   if (zetan) *zetan = z;
     919        7906 :   return a;
     920             : }
     921             : 
     922             : ulong
     923       85834 : Fl_sqrtn(ulong a, long n, ulong p, ulong *zetan)
     924             : {
     925       85834 :   ulong pi = get_Fl_red(p);
     926       85834 :   return Fl_sqrtn_pre(a, n, p, pi, zetan);
     927             : }
     928             : 
     929             : /* Cipolla is better than Tonelli-Shanks when e = v_2(p-1) is "too big".
     930             :  * Otherwise, is a constant times worse; for p = 3 (mod 4), is about 3 times worse,
     931             :  * and in average is about 2 or 2.5 times worse. But try both algorithms for
     932             :  * S(n) = (2^n+3)^2-8 with n = 750, 771, 779, 790, 874, 1176, 1728, 2604, etc.
     933             :  *
     934             :  * If X^2 := t^2 - a  is not a square in F_p (so X is in F_p^2), then
     935             :  *   (t+X)^(p+1) = (t-X)(t+X) = a,   hence  sqrt(a) = (t+X)^((p+1)/2)  in F_p^2.
     936             :  * If (a|p)=1, then sqrt(a) is in F_p.
     937             :  * cf: LNCS 2286, pp 430-434 (2002)  [Gonzalo Tornaria] */
     938             : 
     939             : /* compute y^2, y = y[1] + y[2] X */
     940             : static GEN
     941         449 : sqrt_Cipolla_sqr(void *data, GEN y)
     942             : {
     943         449 :   GEN u = gel(y,1), v = gel(y,2), p = gel(data,2), n = gel(data,3);
     944         449 :   GEN u2 = sqri(u), v2 = sqri(v);
     945         449 :   v = subii(sqri(addii(v,u)), addii(u2,v2));
     946         449 :   u = addii(u2, mulii(v2,n));
     947         449 :   retmkvec2(modii(u,p), modii(v,p));
     948             : }
     949             : /* compute (t+X) y^2 */
     950             : static GEN
     951          23 : sqrt_Cipolla_msqr(void *data, GEN y)
     952             : {
     953          23 :   GEN u = gel(y,1), v = gel(y,2), a = gel(data,1), p = gel(data,2);
     954          23 :   ulong t = gel(data,4)[2];
     955          23 :   GEN d = addii(u, mului(t,v)), d2 = sqri(d);
     956          23 :   GEN b = remii(mulii(a,v), p);
     957          23 :   u = subii(mului(t,d2), mulii(b,addii(u,d)));
     958          23 :   v = subii(d2, mulii(b,v));
     959          23 :   retmkvec2(modii(u,p), modii(v,p));
     960             : }
     961             : /* assume a reduced mod p [ otherwise correct but inefficient ] */
     962             : static GEN
     963           8 : sqrt_Cipolla(GEN a, GEN p)
     964             : {
     965             :   pari_sp av;
     966             :   GEN u, v, n, y, pov2;
     967             :   ulong t;
     968             : 
     969           8 :   if (kronecker(a, p) < 0) return NULL;
     970           8 :   pov2 = shifti(p,-1); /* center to avoid multiplying by huge base*/
     971           8 :   if (cmpii(a,pov2) > 0) a = subii(a,p);
     972           8 :   av = avma;
     973          41 :   for (t=1; ; t++, set_avma(av))
     974             :   {
     975          41 :     n = subsi((long)(t*t), a);
     976          41 :     if (kronecker(n, p) < 0) break;
     977             :   }
     978             : 
     979             :   /* compute (t+X)^((p-1)/2) =: u+vX */
     980           8 :   u = utoipos(t);
     981           8 :   y = gen_pow_fold(mkvec2(u, gen_1), pov2, mkvec4(a,p,n,u),
     982             :                    sqrt_Cipolla_sqr, sqrt_Cipolla_msqr);
     983             :   /* Now u+vX = (t+X)^((p-1)/2); thus
     984             :    *   (u+vX)(t+X) = sqrt(a) + 0 X
     985             :    * Whence,
     986             :    *   sqrt(a) = (u+vt)t - v*a
     987             :    *   0       = (u+vt)
     988             :    * Thus a square root is v*a */
     989             : 
     990           8 :   v = Fp_mul(gel(y, 2), a, p);
     991           8 :   if (cmpii(v,pov2) > 0) v = subii(p,v);
     992           8 :   return v;
     993             : }
     994             : 
     995             : /* Return NULL if p is found to be composite */
     996             : static GEN
     997        3190 : Fp_2gener_all(long e, GEN p)
     998             : {
     999             :   GEN y, m;
    1000             :   long k;
    1001        3190 :   GEN q = shifti(subiu(p,1), -e); /* q = (p-1)/2^oo is odd */
    1002        3190 :   if (e==0 && !equaliu(p,2)) return NULL;
    1003        3190 :   for (k=2; ; k++)
    1004        7811 :   {
    1005       11001 :     long i = krosi(k, p);
    1006       11001 :     if (i >= 0)
    1007             :     {
    1008        7811 :       if (i) continue;
    1009           0 :       return NULL;
    1010             :     }
    1011        3190 :     y = m = Fp_pow(utoi(k), q, p);
    1012       10928 :     for (i=1; i<e; i++)
    1013        7738 :       if (equali1(m = Fp_sqr(m, p))) break;
    1014        3190 :     if (i == e) break; /* success */
    1015             :   }
    1016        3190 :   return y;
    1017             : }
    1018             : 
    1019             : /* Return NULL if p is found to be composite */
    1020             : GEN
    1021        1120 : Fp_2gener(GEN p)
    1022        1120 : { return Fp_2gener_all(vali(subis(p,1)),p); }
    1023             : 
    1024             : /* smallest square root */
    1025             : static GEN
    1026       32996 : choose_sqrt(GEN v, GEN p)
    1027             : {
    1028       32996 :   pari_sp av = avma;
    1029       32996 :   GEN q = subii(p,v);
    1030       32995 :   if (cmpii(v,q) > 0) v = q; else set_avma(av);
    1031       32995 :   return v;
    1032             : }
    1033             : /* Tonelli-Shanks. Assume p is prime and return NULL if (a,p) = -1. */
    1034             : GEN
    1035     3334227 : Fp_sqrt_i(GEN a, GEN y, GEN p)
    1036             : {
    1037     3334227 :   pari_sp av = avma;
    1038             :   long i, k, e;
    1039             :   GEN p1, q, v, w;
    1040             : 
    1041     3334227 :   if (lgefint(p) == 3)
    1042             :   {
    1043     3301126 :     ulong pp = uel(p,2), u = umodiu(a, pp);
    1044     3301157 :     if (!u) return gen_0;
    1045     3005441 :     u = Fl_sqrt(u, pp);
    1046     3005501 :     return (u == ~0UL)? NULL: utoipos(u);
    1047             :   }
    1048             : 
    1049       33101 :   a = modii(a, p); if (!signe(a)) return gc_const(av, gen_0);
    1050       32984 :   p1 = subiu(p,1); e = vali(p1);
    1051       32998 :   if (e <= 2)
    1052             :   { /* direct formulas more efficient */
    1053             :     pari_sp av2;
    1054       26519 :     if (e == 0) pari_err_PRIME("Fp_sqrt [modulus]",p); /* p != 2 */
    1055       26519 :     if (e == 1)
    1056             :     {
    1057       16011 :       q = addiu(shifti(p1,-2),1); /* (p+1) / 4 */
    1058       16003 :       v = Fp_pow(a, q, p);
    1059             :     }
    1060             :     else
    1061             :     { /* Atkin's formula */
    1062       10508 :       GEN i, a2 = shifti(a,1);
    1063       10506 :       if (cmpii(a2,p) >= 0) a2 = subii(a2,p);
    1064       10510 :       q = shifti(p1, -3); /* (p-5)/8 */
    1065       10512 :       v = Fp_pow(a2, q, p);
    1066       10519 :       i = Fp_mul(a2, Fp_sqr(v,p), p); /* i^2 = -1 */
    1067       10519 :       v = Fp_mul(a, Fp_mul(v, subiu(i,1), p), p);
    1068             :     }
    1069       26545 :     av2 = avma;
    1070             :     /* must check equality in case (a/p) = -1 or p not prime */
    1071       26545 :     e = equalii(Fp_sqr(v,p), a); set_avma(av2);
    1072       26546 :     return e? gerepileuptoint(av,choose_sqrt(v,p)): NULL;
    1073             :   }
    1074             :   /* On average, Cipolla is better than Tonelli/Shanks if and only if
    1075             :    * e(e-1) > 8*log2(n)+20, see LNCS 2286 pp 430 [GTL] */
    1076        6479 :   if (e*(e-1) > 20 + 8 * expi(p))
    1077             :   {
    1078           8 :     v = sqrt_Cipolla(a,p); if (!v) return gc_NULL(av);
    1079           8 :     return gerepileuptoint(av,v);
    1080             :   }
    1081        6472 :   if (!y)
    1082             :   {
    1083        2070 :     y = Fp_2gener_all(e, p);
    1084        2070 :     if (!y) pari_err_PRIME("Fp_sqrt [modulus]",p);
    1085             :   }
    1086        6472 :   q = shifti(p1,-e); /* q = (p-1)/2^oo is odd */
    1087        6471 :   p1 = Fp_pow(a, shifti(q,-1), p); /* a ^ (q-1)/2 */
    1088        6472 :   v = Fp_mul(a, p1, p);
    1089        6473 :   w = Fp_mul(v, p1, p);
    1090       15684 :   while (!equali1(w))
    1091             :   { /* a*w = v^2, y primitive 2^e-th root of 1
    1092             :        a square --> w even power of y, hence w^(2^(e-1)) = 1 */
    1093        9212 :     p1 = Fp_sqr(w,p);
    1094       19896 :     for (k=1; !equali1(p1) && k < e; k++) p1 = Fp_sqr(p1,p);
    1095        9204 :     if (k == e) return gc_NULL(av); /* p composite or (a/p) != 1 */
    1096             :     /* w ^ (2^k) = 1 --> w = y ^ (u * 2^(e-k)), u odd */
    1097        9204 :     p1 = y;
    1098       13380 :     for (i=1; i < e-k; i++) p1 = Fp_sqr(p1,p);
    1099        9204 :     y = Fp_sqr(p1, p); e = k;
    1100        9211 :     w = Fp_mul(y, w, p);
    1101        9212 :     v = Fp_mul(v, p1, p);
    1102        9211 :     if (gc_needed(av,1))
    1103             :     {
    1104           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"Fp_sqrt");
    1105           0 :       gerepileall(av,3, &y,&w,&v);
    1106             :     }
    1107             :   }
    1108        6470 :   return gerepileuptoint(av, choose_sqrt(v,p));
    1109             : }
    1110             : 
    1111             : GEN
    1112     3314854 : Fp_sqrt(GEN a, GEN p)
    1113     3314854 : { return Fp_sqrt_i(a, NULL, p); }
    1114             : 
    1115             : /*********************************************************************/
    1116             : /**                        GCD & BEZOUT                             **/
    1117             : /*********************************************************************/
    1118             : 
    1119             : GEN
    1120    41157101 : lcmii(GEN x, GEN y)
    1121             : {
    1122             :   pari_sp av;
    1123             :   GEN a, b;
    1124    41157101 :   if (!signe(x) || !signe(y)) return gen_0;
    1125    41157192 :   av = avma; a = gcdii(x,y);
    1126    41155654 :   if (absequalii(a,y)) { set_avma(av); return absi(x); }
    1127     7805305 :   if (!equali1(a)) y = diviiexact(y,a);
    1128     7805343 :   b = mulii(x,y); setabssign(b); return gerepileuptoint(av, b);
    1129             : }
    1130             : 
    1131             : /* given x in assume 0 < x < N; return u in (Z/NZ)^* such that u x = gcd(x,N) (mod N);
    1132             :  * set *pd = gcd(x,N) */
    1133             : GEN
    1134     6587193 : Fp_invgen(GEN x, GEN N, GEN *pd)
    1135             : {
    1136             :   GEN d, d0, e, v;
    1137     6587193 :   if (lgefint(N) == 3)
    1138             :   {
    1139     6026489 :     ulong dd, NN = N[2], xx = umodiu(x,NN);
    1140     6027123 :     if (!xx) { *pd = N; return gen_0; }
    1141     6027123 :     xx = Fl_invgen(xx, NN, &dd);
    1142     6027718 :     *pd = utoi(dd); return utoi(xx);
    1143             :   }
    1144      560704 :   *pd = d = bezout(x, N, &v, NULL);
    1145      560717 :   if (equali1(d)) return v;
    1146             :   /* vx = gcd(x,N) (mod N), v coprime to N/d but need not be coprime to N */
    1147      460833 :   e = diviiexact(N,d);
    1148      460833 :   d0 = Z_ppo(d, e); /* d = d0 d1, d0 coprime to N/d, rad(d1) | N/d */
    1149      460833 :   if (equali1(d0)) return v;
    1150      289792 :   if (!equalii(d,d0)) e = lcmii(e, diviiexact(d,d0));
    1151      289792 :   return Z_chinese_coprime(v, gen_1, e, d0, mulii(e,d0));
    1152             : }
    1153             : 
    1154             : /*********************************************************************/
    1155             : /**                      CHINESE REMAINDERS                         **/
    1156             : /*********************************************************************/
    1157             : 
    1158             : /* Chinese Remainder Theorem.  x and y must have the same type (integermod,
    1159             :  * polymod, or polynomial/vector/matrix recursively constructed with these
    1160             :  * as coefficients). Creates (with the same type) a z in the same residue
    1161             :  * class as x and the same residue class as y, if it is possible.
    1162             :  *
    1163             :  * We also allow (during recursion) two identical objects even if they are
    1164             :  * not integermod or polymod. For example:
    1165             :  *
    1166             :  * ? x = [1, Mod(5, 11), Mod(X + Mod(2, 7), X^2 + 1)];
    1167             :  * ? y = [1, Mod(7, 17), Mod(X + Mod(0, 3), X^2 + 1)];
    1168             :  * ? chinese(x, y)
    1169             :  * %3 = [1, Mod(16, 187), Mod(X + mod(9, 21), X^2 + 1)] */
    1170             : 
    1171             : static GEN
    1172     2127668 : gen_chinese(GEN x, GEN(*f)(GEN,GEN))
    1173             : {
    1174     2127668 :   GEN z = gassoc_proto(f,x,NULL);
    1175     2127661 :   if (z == gen_1) retmkintmod(gen_0,gen_1);
    1176     2127626 :   return z;
    1177             : }
    1178             : 
    1179             : /* x t_INTMOD, y t_POLMOD; promote x to t_POLMOD mod Pol(x.mod) then
    1180             :  * call chinese: makes Mod(0,1) a better "neutral" element */
    1181             : static GEN
    1182          21 : chinese_intpol(GEN x,GEN y)
    1183             : {
    1184          21 :   pari_sp av = avma;
    1185          21 :   GEN z = mkpolmod(gel(x,2), scalarpol_shallow(gel(x,1), varn(gel(y,1))));
    1186          21 :   return gerepileupto(av, chinese(z, y));
    1187             : }
    1188             : 
    1189             : GEN
    1190        2275 : chinese1(GEN x) { return gen_chinese(x,chinese); }
    1191             : 
    1192             : GEN
    1193        5327 : chinese(GEN x, GEN y)
    1194             : {
    1195             :   pari_sp av;
    1196        5327 :   long tx = typ(x), ty;
    1197             :   GEN z,p1,p2,d,u,v;
    1198             : 
    1199        5327 :   if (!y) return chinese1(x);
    1200        5278 :   if (gidentical(x,y)) return gcopy(x);
    1201        5271 :   ty = typ(y);
    1202        5271 :   if (tx == ty) switch(tx)
    1203             :   {
    1204        3759 :     case t_POLMOD:
    1205             :     {
    1206        3759 :       GEN A = gel(x,1), B = gel(y,1);
    1207        3759 :       GEN a = gel(x,2), b = gel(y,2);
    1208        3759 :       if (varn(A)!=varn(B)) pari_err_VAR("chinese",A,B);
    1209        3759 :       if (RgX_equal(A,B)) retmkpolmod(chinese(a,b), gcopy(A)); /*same modulus*/
    1210        3759 :       av = avma;
    1211        3759 :       d = RgX_extgcd(A,B,&u,&v);
    1212        3759 :       p2 = gsub(b, a);
    1213        3759 :       if (!gequal0(gmod(p2, d))) break;
    1214        3759 :       p1 = gdiv(A,d);
    1215        3759 :       p2 = gadd(a, gmul(gmul(u,p1), p2));
    1216             : 
    1217        3759 :       z = cgetg(3, t_POLMOD);
    1218        3759 :       gel(z,1) = gmul(p1,B);
    1219        3759 :       gel(z,2) = gmod(p2,gel(z,1));
    1220        3759 :       return gerepileupto(av, z);
    1221             :     }
    1222        1477 :     case t_INTMOD:
    1223             :     {
    1224        1477 :       GEN A = gel(x,1), B = gel(y,1);
    1225        1477 :       GEN a = gel(x,2), b = gel(y,2), c, d, C, U;
    1226        1477 :       z = cgetg(3,t_INTMOD);
    1227        1477 :       Z_chinese_pre(A, B, &C, &U, &d);
    1228        1477 :       c = Z_chinese_post(a, b, C, U, d);
    1229        1477 :       if (!c) pari_err_OP("chinese", x,y);
    1230        1477 :       set_avma((pari_sp)z);
    1231        1477 :       gel(z,1) = icopy(C);
    1232        1477 :       gel(z,2) = icopy(c); return z;
    1233             :     }
    1234           7 :     case t_POL:
    1235             :     {
    1236           7 :       long i, lx = lg(x), ly = lg(y);
    1237           7 :       if (varn(x) != varn(y)) break;
    1238           7 :       if (lx < ly) { swap(x,y); lswap(lx,ly); }
    1239           7 :       z = cgetg(lx, t_POL); z[1] = x[1];
    1240          21 :       for (i=2; i<ly; i++) gel(z,i) = chinese(gel(x,i),gel(y,i));
    1241          14 :       for (   ; i<lx; i++) gel(z,i) = gcopy(gel(x,i));
    1242           7 :       return z;
    1243             :     }
    1244             : 
    1245           7 :     case t_VEC: case t_COL: case t_MAT:
    1246             :     {
    1247             :       long i, lx;
    1248           7 :       z = cgetg_copy(x, &lx); if (lx!=lg(y)) break;
    1249          21 :       for (i=1; i<lx; i++) gel(z,i) = chinese(gel(x,i),gel(y,i));
    1250           7 :       return z;
    1251             :     }
    1252             :   }
    1253          21 :   if (tx == t_POLMOD && ty == t_INTMOD) return chinese_intpol(y,x);
    1254           7 :   if (ty == t_POLMOD && tx == t_INTMOD) return chinese_intpol(x,y);
    1255           0 :   pari_err_OP("chinese",x,y);
    1256             :   return NULL; /* LCOV_EXCL_LINE */
    1257             : }
    1258             : 
    1259             : /* init chinese(Mod(.,A), Mod(.,B)) */
    1260             : void
    1261      266158 : Z_chinese_pre(GEN A, GEN B, GEN *pC, GEN *pU, GEN *pd)
    1262             : {
    1263      266158 :   GEN u, d = bezout(A,B,&u,NULL); /* U = u(A/d), u(A/d) + v(B/d) = 1 */
    1264      266162 :   GEN t = diviiexact(A,d);
    1265      266153 :   *pU = mulii(u, t);
    1266      266146 :   *pC = mulii(t, B);
    1267      266145 :   if (pd) *pd = d;
    1268      266145 : }
    1269             : /* Assume C = lcm(A, B), U = 0 mod (A/d), U = 1 mod (B/d), a = b mod d,
    1270             :  * where d = gcd(A,B) or NULL, return x = a (mod A), b (mod B).
    1271             :  * If d not NULL, check whether a = b mod d. */
    1272             : GEN
    1273     2730476 : Z_chinese_post(GEN a, GEN b, GEN C, GEN U, GEN d)
    1274             : {
    1275             :   GEN b_a;
    1276     2730476 :   if (!signe(a))
    1277             :   {
    1278      785220 :     if (d && !dvdii(b, d)) return NULL;
    1279      785220 :     return Fp_mul(b, U, C);
    1280             :   }
    1281     1945256 :   b_a = subii(b,a);
    1282     1945256 :   if (d && !dvdii(b_a, d)) return NULL;
    1283     1945256 :   return modii(addii(a, mulii(U, b_a)), C);
    1284             : }
    1285             : static ulong
    1286     2114536 : u_chinese_post(ulong a, ulong b, ulong C, ulong U)
    1287             : {
    1288     2114536 :   if (!a) return Fl_mul(b, U, C);
    1289     2114536 :   return Fl_add(a, Fl_mul(U, Fl_sub(b,a,C), C), C);
    1290             : }
    1291             : 
    1292             : GEN
    1293        2142 : Z_chinese(GEN a, GEN b, GEN A, GEN B)
    1294             : {
    1295        2142 :   pari_sp av = avma;
    1296        2142 :   GEN C, U; Z_chinese_pre(A, B, &C, &U, NULL);
    1297        2142 :   return gerepileuptoint(av, Z_chinese_post(a,b, C, U, NULL));
    1298             : }
    1299             : GEN
    1300      262483 : Z_chinese_all(GEN a, GEN b, GEN A, GEN B, GEN *pC)
    1301             : {
    1302      262483 :   GEN U; Z_chinese_pre(A, B, pC, &U, NULL);
    1303      262469 :   return Z_chinese_post(a,b, *pC, U, NULL);
    1304             : }
    1305             : 
    1306             : /* return lift(chinese(a mod A, b mod B))
    1307             :  * assume(A,B)=1, a,b,A,B integers and C = A*B */
    1308             : GEN
    1309      291052 : Z_chinese_coprime(GEN a, GEN b, GEN A, GEN B, GEN C)
    1310             : {
    1311      291052 :   pari_sp av = avma;
    1312      291052 :   GEN U = mulii(Fp_inv(A,B), A);
    1313      291052 :   return gerepileuptoint(av, Z_chinese_post(a,b,C,U, NULL));
    1314             : }
    1315             : ulong
    1316     2114521 : u_chinese_coprime(ulong a, ulong b, ulong A, ulong B, ulong C)
    1317     2114521 : { return u_chinese_post(a,b,C, A * Fl_inv(A % B,B)); }
    1318             : 
    1319             : /* chinese1 for coprime moduli in Z */
    1320             : static GEN
    1321     2173011 : chinese1_coprime_Z_aux(GEN x, GEN y)
    1322             : {
    1323     2173011 :   GEN z = cgetg(3, t_INTMOD);
    1324     2173011 :   GEN A = gel(x,1), a = gel(x, 2);
    1325     2173011 :   GEN B = gel(y,1), b = gel(y, 2), C = mulii(A,B);
    1326     2173011 :   pari_sp av = avma;
    1327     2173011 :   GEN U = mulii(Fp_inv(A,B), A);
    1328     2173011 :   gel(z,2) = gerepileuptoint(av, Z_chinese_post(a,b,C,U, NULL));
    1329     2173011 :   gel(z,1) = C; return z;
    1330             : }
    1331             : GEN
    1332     2125393 : chinese1_coprime_Z(GEN x) {return gen_chinese(x,chinese1_coprime_Z_aux);}
    1333             : 
    1334             : /*********************************************************************/
    1335             : /**                    MODULAR EXPONENTIATION                       **/
    1336             : /*********************************************************************/
    1337             : /* xa ZV or nv */
    1338             : GEN
    1339     1998471 : ZV_producttree(GEN xa)
    1340             : {
    1341     1998471 :   long n = lg(xa)-1;
    1342     1998471 :   long m = n==1 ? 1: expu(n-1)+1;
    1343     1998470 :   GEN T = cgetg(m+1, t_VEC), t;
    1344             :   long i, j, k;
    1345     1998470 :   t = cgetg(((n+1)>>1)+1, t_VEC);
    1346     1998461 :   if (typ(xa)==t_VECSMALL)
    1347             :   {
    1348     2581316 :     for (j=1, k=1; k<n; j++, k+=2)
    1349     1703417 :       gel(t, j) = muluu(xa[k], xa[k+1]);
    1350      877899 :     if (k==n) gel(t, j) = utoi(xa[k]);
    1351             :   } else {
    1352     2311315 :     for (j=1, k=1; k<n; j++, k+=2)
    1353     1190781 :       gel(t, j) = mulii(gel(xa,k), gel(xa,k+1));
    1354     1120534 :     if (k==n) gel(t, j) = icopy(gel(xa,k));
    1355             :   }
    1356     1998437 :   gel(T,1) = t;
    1357     3269065 :   for (i=2; i<=m; i++)
    1358             :   {
    1359     1270618 :     GEN u = gel(T, i-1);
    1360     1270618 :     long n = lg(u)-1;
    1361     1270618 :     t = cgetg(((n+1)>>1)+1, t_VEC);
    1362     2830025 :     for (j=1, k=1; k<n; j++, k+=2)
    1363     1559397 :       gel(t, j) = mulii(gel(u, k), gel(u, k+1));
    1364     1270628 :     if (k==n) gel(t, j) = gel(u, k);
    1365     1270628 :     gel(T, i) = t;
    1366             :   }
    1367     1998447 :   return T;
    1368             : }
    1369             : 
    1370             : /* return [A mod P[i], i=1..#P], T = ZV_producttree(P) */
    1371             : GEN
    1372    53043330 : Z_ZV_mod_tree(GEN A, GEN P, GEN T)
    1373             : {
    1374             :   long i,j,k;
    1375    53043330 :   long m = lg(T)-1, n = lg(P)-1;
    1376             :   GEN t;
    1377    53043330 :   GEN Tp = cgetg(m+1, t_VEC);
    1378    52980002 :   gel(Tp, m) = mkvec(modii(A, gmael(T,m,1)));
    1379   111639528 :   for (i=m-1; i>=1; i--)
    1380             :   {
    1381    58753897 :     GEN u = gel(T, i);
    1382    58753897 :     GEN v = gel(Tp, i+1);
    1383    58753897 :     long n = lg(u)-1;
    1384    58753897 :     t = cgetg(n+1, t_VEC);
    1385   146307430 :     for (j=1, k=1; k<n; j++, k+=2)
    1386             :     {
    1387    87689976 :       gel(t, k)   = modii(gel(v, j), gel(u, k));
    1388    87720368 :       gel(t, k+1) = modii(gel(v, j), gel(u, k+1));
    1389             :     }
    1390    58617454 :     if (k==n) gel(t, k) = gel(v, j);
    1391    58617454 :     gel(Tp, i) = t;
    1392             :   }
    1393             :   {
    1394    52885631 :     GEN u = gel(T, i+1);
    1395    52885631 :     GEN v = gel(Tp, i+1);
    1396    52885631 :     long l = lg(u)-1;
    1397    52885631 :     if (typ(P)==t_VECSMALL)
    1398             :     {
    1399    50892637 :       GEN R = cgetg(n+1, t_VECSMALL);
    1400   188783109 :       for (j=1, k=1; j<=l; j++, k+=2)
    1401             :       {
    1402   137581366 :         uel(R,k) = umodiu(gel(v, j), P[k]);
    1403   137851467 :         if (k < n)
    1404   110527042 :           uel(R,k+1) = umodiu(gel(v, j), P[k+1]);
    1405             :       }
    1406    51201743 :       return R;
    1407             :     }
    1408             :     else
    1409             :     {
    1410     1992994 :       GEN R = cgetg(n+1, t_VEC);
    1411     5552869 :       for (j=1, k=1; j<=l; j++, k+=2)
    1412             :       {
    1413     3554706 :         gel(R,k) = modii(gel(v, j), gel(P,k));
    1414     3554678 :         if (k < n)
    1415     2891105 :           gel(R,k+1) = modii(gel(v, j), gel(P,k+1));
    1416             :       }
    1417     1998163 :       return R;
    1418             :     }
    1419             :   }
    1420             : }
    1421             : 
    1422             : /* T = ZV_producttree(P), R = ZV_chinesetree(P,T) */
    1423             : GEN
    1424    34778835 : ZV_chinese_tree(GEN A, GEN P, GEN T, GEN R)
    1425             : {
    1426    34778835 :   long m = lg(T)-1, n = lg(A)-1;
    1427             :   long i,j,k;
    1428    34778835 :   GEN Tp = cgetg(m+1, t_VEC);
    1429    34751463 :   GEN M = gel(T, 1);
    1430    34751463 :   GEN t = cgetg(lg(M), t_VEC);
    1431    34709715 :   if (typ(P)==t_VECSMALL)
    1432             :   {
    1433    80887346 :     for (j=1, k=1; k<n; j++, k+=2)
    1434             :     {
    1435    60411866 :       pari_sp av = avma;
    1436    60411866 :       GEN a = mului(A[k], gel(R,k)), b = mului(A[k+1], gel(R,k+1));
    1437    60191831 :       GEN tj = modii(addii(mului(P[k],b), mului(P[k+1],a)), gel(M,j));
    1438    60351224 :       gel(t, j) = gerepileuptoint(av, tj);
    1439             :     }
    1440    20475480 :     if (k==n) gel(t, j) = modii(mului(A[k], gel(R,k)), gel(M, j));
    1441             :   } else
    1442             :   {
    1443    30318165 :     for (j=1, k=1; k<n; j++, k+=2)
    1444             :     {
    1445    16047924 :       pari_sp av = avma;
    1446    16047924 :       GEN a = mulii(gel(A,k), gel(R,k)), b = mulii(gel(A,k+1), gel(R,k+1));
    1447    16059331 :       GEN tj = modii(addii(mulii(gel(P,k),b), mulii(gel(P,k+1),a)), gel(M,j));
    1448    16085307 :       gel(t, j) = gerepileuptoint(av, tj);
    1449             :     }
    1450    14270241 :     if (k==n) gel(t, j) = modii(mulii(gel(A,k), gel(R,k)), gel(M, j));
    1451             :   }
    1452    34736173 :   gel(Tp, 1) = t;
    1453    66811088 :   for (i=2; i<=m; i++)
    1454             :   {
    1455    32060465 :     GEN u = gel(T, i-1), M = gel(T, i);
    1456    32060465 :     GEN t = cgetg(lg(M), t_VEC);
    1457    32092324 :     GEN v = gel(Tp, i-1);
    1458    32092324 :     long n = lg(v)-1;
    1459    87492007 :     for (j=1, k=1; k<n; j++, k+=2)
    1460             :     {
    1461    55417092 :       pari_sp av = avma;
    1462    55357501 :       gel(t, j) = gerepileuptoint(av, modii(addii(mulii(gel(u, k), gel(v, k+1)),
    1463    55417092 :             mulii(gel(u, k+1), gel(v, k))), gel(M, j)));
    1464             :     }
    1465    32074915 :     if (k==n) gel(t, j) = gel(v, k);
    1466    32074915 :     gel(Tp, i) = t;
    1467             :   }
    1468    34750623 :   return gmael(Tp,m,1);
    1469             : }
    1470             : 
    1471             : static GEN
    1472      931544 : ncV_polint_center_tree(GEN vA, GEN P, GEN T, GEN R, GEN m2)
    1473             : {
    1474      931544 :   long i, l = lg(gel(vA,1)), n = lg(P);
    1475      931544 :   GEN mod = gmael(T, lg(T)-1, 1), V = cgetg(l, t_COL);
    1476    30520609 :   for (i=1; i < l; i++)
    1477             :   {
    1478    29589191 :     pari_sp av = avma;
    1479    29589191 :     GEN c, A = cgetg(n, typ(P));
    1480             :     long j;
    1481   181956313 :     for (j=1; j < n; j++) A[j] = mael(vA,j,i);
    1482    29553802 :     c = Fp_center(ZV_chinese_tree(A, P, T, R), mod, m2);
    1483    29588021 :     gel(V,i) = gerepileuptoint(av, c);
    1484             :   }
    1485      931418 :   return V;
    1486             : }
    1487             : 
    1488             : static GEN
    1489      392441 : nxV_polint_center_tree(GEN vA, GEN P, GEN T, GEN R, GEN m2)
    1490             : {
    1491      392441 :   long i, j, l, n = lg(P);
    1492      392441 :   GEN mod = gmael(T, lg(T)-1, 1), V, w;
    1493      392441 :   w = cgetg(n, t_VECSMALL);
    1494     1434443 :   for(j=1; j<n; j++) w[j] = lg(gel(vA,j));
    1495      392428 :   l = vecsmall_max(w);
    1496      392433 :   V = cgetg(l, t_POL);
    1497      392416 :   V[1] = mael(vA,1,1);
    1498     2591519 :   for (i=2; i < l; i++)
    1499             :   {
    1500     2199086 :     pari_sp av = avma;
    1501     2199086 :     GEN c, A = cgetg(n, typ(P));
    1502     2198803 :     if (typ(P)==t_VECSMALL)
    1503     4669270 :       for (j=1; j < n; j++) A[j] = i < w[j] ? mael(vA,j,i): 0;
    1504             :     else
    1505     3739003 :       for (j=1; j < n; j++) gel(A,j) = i < w[j] ? gmael(vA,j,i): gen_0;
    1506     2198803 :     c = Fp_center(ZV_chinese_tree(A, P, T, R), mod, m2);
    1507     2199190 :     gel(V,i) = gerepileuptoint(av, c);
    1508             :   }
    1509      392433 :   return ZX_renormalize(V, l);
    1510             : }
    1511             : 
    1512             : static GEN
    1513        4614 : nxCV_polint_center_tree(GEN vA, GEN P, GEN T, GEN R, GEN m2)
    1514             : {
    1515        4614 :   long i, j, l = lg(gel(vA,1)), n = lg(P);
    1516        4614 :   GEN A = cgetg(n, t_VEC);
    1517        4614 :   GEN V = cgetg(l, t_COL);
    1518       90887 :   for (i=1; i < l; i++)
    1519             :   {
    1520      335054 :     for (j=1; j < n; j++) gel(A,j) = gmael(vA,j,i);
    1521       86273 :     gel(V,i) = nxV_polint_center_tree(A, P, T, R, m2);
    1522             :   }
    1523        4614 :   return V;
    1524             : }
    1525             : 
    1526             : static GEN
    1527      114633 : polint_chinese(GEN worker, GEN mA, GEN P)
    1528             : {
    1529      114633 :   long cnt, pending, n, i, j, l = lg(gel(mA,1));
    1530             :   struct pari_mt pt;
    1531             :   GEN done, va, M, A;
    1532             :   pari_timer ti;
    1533             : 
    1534      114633 :   if (l == 1) return cgetg(1, t_MAT);
    1535       85581 :   cnt = pending = 0;
    1536       85581 :   n = lg(P);
    1537       85581 :   A = cgetg(n, t_VEC);
    1538       85581 :   va = mkvec(A);
    1539       85581 :   M = cgetg(l, t_MAT);
    1540       85581 :   if (DEBUGLEVEL>4) timer_start(&ti);
    1541       85581 :   if (DEBUGLEVEL>5) err_printf("Start parallel Chinese remainder: ");
    1542       85581 :   mt_queue_start_lim(&pt, worker, l-1);
    1543      671529 :   for (i=1; i<l || pending; i++)
    1544             :   {
    1545             :     long workid;
    1546     2421256 :     for(j=1; j < n; j++) gel(A,j) = gmael(mA,j,i);
    1547      585948 :     mt_queue_submit(&pt, i, i<l? va: NULL);
    1548      585948 :     done = mt_queue_get(&pt, &workid, &pending);
    1549      585948 :     if (done)
    1550             :     {
    1551      552208 :       gel(M,workid) = done;
    1552      552208 :       if (DEBUGLEVEL>5) err_printf("%ld%% ",(++cnt)*100/(l-1));
    1553             :     }
    1554             :   }
    1555       85581 :   if (DEBUGLEVEL>5) err_printf("\n");
    1556       85581 :   if (DEBUGLEVEL>4) timer_printf(&ti, "nmV_chinese_center");
    1557       85581 :   mt_queue_end(&pt);
    1558       85581 :   return M;
    1559             : }
    1560             : 
    1561             : GEN
    1562         840 : nxMV_polint_center_tree_worker(GEN vA, GEN T, GEN R, GEN P, GEN m2)
    1563             : {
    1564         840 :   return nxCV_polint_center_tree(vA, P, T, R, m2);
    1565             : }
    1566             : 
    1567             : static GEN
    1568         430 : nxMV_polint_center_tree_seq(GEN vA, GEN P, GEN T, GEN R, GEN m2)
    1569             : {
    1570         430 :   long i, j, l = lg(gel(vA,1)), n = lg(P);
    1571         430 :   GEN A = cgetg(n, t_VEC);
    1572         430 :   GEN V = cgetg(l, t_MAT);
    1573        4204 :   for (i=1; i < l; i++)
    1574             :   {
    1575       15299 :     for (j=1; j < n; j++) gel(A,j) = gmael(vA,j,i);
    1576        3774 :     gel(V,i) = nxCV_polint_center_tree(A, P, T, R, m2);
    1577             :   }
    1578         430 :   return V;
    1579             : }
    1580             : 
    1581             : static GEN
    1582          90 : nxMV_polint_center_tree(GEN mA, GEN P, GEN T, GEN R, GEN m2)
    1583             : {
    1584          90 :   GEN worker = snm_closure(is_entry("_nxMV_polint_worker"), mkvec4(T, R, P, m2));
    1585          90 :   return polint_chinese(worker, mA, P);
    1586             : }
    1587             : 
    1588             : static GEN
    1589       52973 : nmV_polint_center_tree_seq(GEN vA, GEN P, GEN T, GEN R, GEN m2)
    1590             : {
    1591       52973 :   long i, j, l = lg(gel(vA,1)), n = lg(P);
    1592       52973 :   GEN A = cgetg(n, t_VEC);
    1593       52973 :   GEN V = cgetg(l, t_MAT);
    1594      418874 :   for (i=1; i < l; i++)
    1595             :   {
    1596     2132700 :     for (j=1; j < n; j++) gel(A,j) = gmael(vA,j,i);
    1597      365898 :     gel(V,i) = ncV_polint_center_tree(A, P, T, R, m2);
    1598             :   }
    1599       52976 :   return V;
    1600             : }
    1601             : 
    1602             : GEN
    1603      551268 : nmV_polint_center_tree_worker(GEN vA, GEN T, GEN R, GEN P, GEN m2)
    1604             : {
    1605      551268 :   return ncV_polint_center_tree(vA, P, T, R, m2);
    1606             : }
    1607             : 
    1608             : static GEN
    1609      114543 : nmV_polint_center_tree(GEN mA, GEN P, GEN T, GEN R, GEN m2)
    1610             : {
    1611      114543 :   GEN worker = snm_closure(is_entry("_polint_worker"), mkvec4(T, R, P, m2));
    1612      114543 :   return polint_chinese(worker, mA, P);
    1613             : }
    1614             : 
    1615             : /* return [A mod P[i], i=1..#P] */
    1616             : GEN
    1617           0 : Z_ZV_mod(GEN A, GEN P)
    1618             : {
    1619           0 :   pari_sp av = avma;
    1620           0 :   return gerepilecopy(av, Z_ZV_mod_tree(A, P, ZV_producttree(P)));
    1621             : }
    1622             : /* P a t_VECSMALL */
    1623             : GEN
    1624           0 : Z_nv_mod(GEN A, GEN P)
    1625             : {
    1626           0 :   pari_sp av = avma;
    1627           0 :   return gerepileuptoleaf(av, Z_ZV_mod_tree(A, P, ZV_producttree(P)));
    1628             : }
    1629             : /* B a ZX, T = ZV_producttree(P) */
    1630             : GEN
    1631     1239547 : ZX_nv_mod_tree(GEN B, GEN A, GEN T)
    1632             : {
    1633             :   pari_sp av;
    1634     1239547 :   long i, j, l = lg(B), n = lg(A)-1;
    1635     1239547 :   GEN V = cgetg(n+1, t_VEC);
    1636     6113346 :   for (j=1; j <= n; j++)
    1637             :   {
    1638     4873995 :     gel(V, j) = cgetg(l, t_VECSMALL);
    1639     4873813 :     mael(V, j, 1) = B[1]&VARNBITS;
    1640             :   }
    1641     1239351 :   av = avma;
    1642    12013888 :   for (i=2; i < l; i++)
    1643             :   {
    1644    10775780 :     GEN v = Z_ZV_mod_tree(gel(B, i), A, T);
    1645    74818998 :     for (j=1; j <= n; j++)
    1646    64052394 :       mael(V, j, i) = v[j];
    1647    10766604 :     set_avma(av);
    1648             :   }
    1649     6112621 :   for (j=1; j <= n; j++)
    1650     4874435 :     (void) Flx_renormalize(gel(V, j), l);
    1651     1238186 :   return V;
    1652             : }
    1653             : 
    1654             : static GEN
    1655      159883 : to_ZX(GEN a, long v) { return typ(a)==t_INT? scalarpol(a,v): a; }
    1656             : 
    1657             : GEN
    1658       10277 : ZXX_nv_mod_tree(GEN P, GEN xa, GEN T, long w)
    1659             : {
    1660       10277 :   pari_sp av = avma;
    1661       10277 :   long i, j, l = lg(P), n = lg(xa)-1;
    1662       10277 :   GEN V = cgetg(n+1, t_VEC);
    1663       39626 :   for (j=1; j <= n; j++)
    1664             :   {
    1665       29349 :     gel(V, j) = cgetg(l, t_POL);
    1666       29349 :     mael(V, j, 1) = P[1]&VARNBITS;
    1667             :   }
    1668       88474 :   for (i=2; i < l; i++)
    1669             :   {
    1670       78198 :     GEN v = ZX_nv_mod_tree(to_ZX(gel(P, i), w), xa, T);
    1671      321139 :     for (j=1; j <= n; j++)
    1672      242942 :       gmael(V, j, i) = gel(v,j);
    1673             :   }
    1674       39625 :   for (j=1; j <= n; j++)
    1675       29349 :     (void) FlxX_renormalize(gel(V, j), l);
    1676       10276 :   return gerepilecopy(av, V);
    1677             : }
    1678             : 
    1679             : GEN
    1680        4047 : ZXC_nv_mod_tree(GEN C, GEN xa, GEN T, long w)
    1681             : {
    1682        4047 :   pari_sp av = avma;
    1683        4047 :   long i, j, l = lg(C), n = lg(xa)-1;
    1684        4047 :   GEN V = cgetg(n+1, t_VEC);
    1685       16915 :   for (j = 1; j <= n; j++)
    1686       12870 :     gel(V, j) = cgetg(l, t_COL);
    1687       85716 :   for (i = 1; i < l; i++)
    1688             :   {
    1689       81683 :     GEN v = ZX_nv_mod_tree(to_ZX(gel(C, i), w), xa, T);
    1690      358530 :     for (j = 1; j <= n; j++)
    1691      276859 :       gmael(V, j, i) = gel(v,j);
    1692             :   }
    1693        4033 :   return gerepilecopy(av, V);
    1694             : }
    1695             : 
    1696             : GEN
    1697         430 : ZXM_nv_mod_tree(GEN M, GEN xa, GEN T, long w)
    1698             : {
    1699         430 :   pari_sp av = avma;
    1700         430 :   long i, j, l = lg(M), n = lg(xa)-1;
    1701         430 :   GEN V = cgetg(n+1, t_VEC);
    1702        2083 :   for (j=1; j <= n; j++)
    1703        1653 :     gel(V, j) = cgetg(l, t_MAT);
    1704        4204 :   for (i=1; i < l; i++)
    1705             :   {
    1706        3774 :     GEN v = ZXC_nv_mod_tree(gel(M, i), xa, T, w);
    1707       15299 :     for (j=1; j <= n; j++)
    1708       11525 :       gmael(V, j, i) = gel(v,j);
    1709             :   }
    1710         430 :   return gerepilecopy(av, V);
    1711             : }
    1712             : 
    1713             : GEN
    1714      936200 : ZV_nv_mod_tree(GEN B, GEN A, GEN T)
    1715             : {
    1716             :   pari_sp av;
    1717      936200 :   long i, j, l = lg(B), n = lg(A)-1;
    1718      936200 :   GEN V = cgetg(n+1, t_VEC);
    1719     4637043 :   for (j=1; j <= n; j++) gel(V, j) = cgetg(l, t_VECSMALL);
    1720      936127 :   av = avma;
    1721    41199659 :   for (i=1; i < l; i++)
    1722             :   {
    1723    40267801 :     GEN v = Z_ZV_mod_tree(gel(B, i), A, T);
    1724   226384133 :     for (j=1; j <= n; j++) mael(V, j, i) = v[j];
    1725    40194472 :     set_avma(av);
    1726             :   }
    1727      931858 :   return V;
    1728             : }
    1729             : 
    1730             : GEN
    1731       80586 : ZM_nv_mod_tree(GEN M, GEN xa, GEN T)
    1732             : {
    1733       80586 :   pari_sp av = avma;
    1734       80586 :   long i, j, l = lg(M), n = lg(xa)-1;
    1735       80586 :   GEN V = cgetg(n+1, t_VEC);
    1736      315207 :   for (j=1; j <= n; j++) gel(V, j) = cgetg(l, t_MAT);
    1737     1016552 :   for (i=1; i < l; i++)
    1738             :   {
    1739      935974 :     GEN v = ZV_nv_mod_tree(gel(M, i), xa, T);
    1740     4637166 :     for (j=1; j <= n; j++) gmael(V, j, i) = gel(v,j);
    1741             :   }
    1742       80578 :   return gerepilecopy(av, V);
    1743             : }
    1744             : 
    1745             : static GEN
    1746     1994331 : ZV_sqr(GEN z)
    1747             : {
    1748     1994331 :   long i,l = lg(z);
    1749     1994331 :   GEN x = cgetg(l, t_VEC);
    1750     1994353 :   if (typ(z)==t_VECSMALL)
    1751     4591041 :     for (i=1; i<l; i++) gel(x,i) = sqru(z[i]);
    1752             :   else
    1753     3823221 :     for (i=1; i<l; i++) gel(x,i) = sqri(gel(z,i));
    1754     1994273 :   return x;
    1755             : }
    1756             : 
    1757             : static GEN
    1758    10538062 : ZT_sqr(GEN x)
    1759             : {
    1760    10538062 :   if (typ(x) == t_INT) return sqri(x);
    1761    13796562 :   pari_APPLY_type(t_VEC, ZT_sqr(gel(x,i)))
    1762             : }
    1763             : 
    1764             : static GEN
    1765     1994323 : ZV_invdivexact(GEN y, GEN x)
    1766             : {
    1767     1994323 :   long i, l = lg(y);
    1768     1994323 :   GEN z = cgetg(l,t_VEC);
    1769     1994344 :   if (typ(x)==t_VECSMALL)
    1770     4590603 :     for (i=1; i<l; i++)
    1771             :     {
    1772     3713014 :       pari_sp av = avma;
    1773     3713014 :       ulong a = Fl_inv(umodiu(diviuexact(gel(y,i),x[i]), x[i]), x[i]);
    1774     3713295 :       set_avma(av); gel(z,i) = utoi(a);
    1775             :     }
    1776             :   else
    1777     3823317 :     for (i=1; i<l; i++)
    1778     2706588 :       gel(z,i) = Fp_inv(diviiexact(gel(y,i), gel(x,i)), gel(x,i));
    1779     1994318 :   return z;
    1780             : }
    1781             : 
    1782             : /* P t_VECSMALL or t_VEC of t_INT  */
    1783             : GEN
    1784     1994305 : ZV_chinesetree(GEN P, GEN T)
    1785             : {
    1786     1994305 :   GEN T2 = ZT_sqr(T), P2 = ZV_sqr(P);
    1787     1994310 :   GEN mod = gmael(T,lg(T)-1,1);
    1788     1994310 :   return ZV_invdivexact(Z_ZV_mod_tree(mod, P2, T2), P);
    1789             : }
    1790             : 
    1791             : static GEN
    1792      485179 : gc_chinese(pari_sp av, GEN T, GEN a, GEN *pt_mod)
    1793             : {
    1794      485179 :   if (!pt_mod)
    1795       11159 :     return gerepileupto(av, a);
    1796             :   else
    1797             :   {
    1798      474020 :     GEN mod = gmael(T, lg(T)-1, 1);
    1799      474020 :     gerepileall(av, 2, &a, &mod);
    1800      474021 :     *pt_mod = mod;
    1801      474021 :     return a;
    1802             :   }
    1803             : }
    1804             : 
    1805             : GEN
    1806      160311 : ZV_chinese_center(GEN A, GEN P, GEN *pt_mod)
    1807             : {
    1808      160311 :   pari_sp av = avma;
    1809      160311 :   GEN T = ZV_producttree(P);
    1810      160311 :   GEN R = ZV_chinesetree(P, T);
    1811      160311 :   GEN a = ZV_chinese_tree(A, P, T, R);
    1812      160311 :   GEN mod = gmael(T, lg(T)-1, 1);
    1813      160311 :   GEN ca = Fp_center(a, mod, shifti(mod,-1));
    1814      160311 :   return gc_chinese(av, T, ca, pt_mod);
    1815             : }
    1816             : 
    1817             : GEN
    1818        5088 : ZV_chinese(GEN A, GEN P, GEN *pt_mod)
    1819             : {
    1820        5088 :   pari_sp av = avma;
    1821        5088 :   GEN T = ZV_producttree(P);
    1822        5088 :   GEN R = ZV_chinesetree(P, T);
    1823        5088 :   GEN a = ZV_chinese_tree(A, P, T, R);
    1824        5088 :   return gc_chinese(av, T, a, pt_mod);
    1825             : }
    1826             : 
    1827             : GEN
    1828      110134 : nxV_chinese_center_tree(GEN A, GEN P, GEN T, GEN R)
    1829             : {
    1830      110134 :   pari_sp av = avma;
    1831      110134 :   GEN m2 = shifti(gmael(T, lg(T)-1, 1), -1);
    1832      110134 :   GEN a = nxV_polint_center_tree(A, P, T, R, m2);
    1833      110135 :   return gerepileupto(av, a);
    1834             : }
    1835             : 
    1836             : GEN
    1837      196024 : nxV_chinese_center(GEN A, GEN P, GEN *pt_mod)
    1838             : {
    1839      196024 :   pari_sp av = avma;
    1840      196024 :   GEN T = ZV_producttree(P);
    1841      196023 :   GEN R = ZV_chinesetree(P, T);
    1842      196022 :   GEN m2 = shifti(gmael(T, lg(T)-1, 1), -1);
    1843      196022 :   GEN a = nxV_polint_center_tree(A, P, T, R, m2);
    1844      196024 :   return gc_chinese(av, T, a, pt_mod);
    1845             : }
    1846             : 
    1847             : GEN
    1848        9123 : ncV_chinese_center(GEN A, GEN P, GEN *pt_mod)
    1849             : {
    1850        9123 :   pari_sp av = avma;
    1851        9123 :   GEN T = ZV_producttree(P);
    1852        9123 :   GEN R = ZV_chinesetree(P, T);
    1853        9123 :   GEN m2 = shifti(gmael(T, lg(T)-1, 1), -1);
    1854        9123 :   GEN a = ncV_polint_center_tree(A, P, T, R, m2);
    1855        9123 :   return gc_chinese(av, T, a, pt_mod);
    1856             : }
    1857             : 
    1858             : GEN
    1859        5262 : ncV_chinese_center_tree(GEN A, GEN P, GEN T, GEN R)
    1860             : {
    1861        5262 :   pari_sp av = avma;
    1862        5262 :   GEN m2 = shifti(gmael(T, lg(T)-1, 1), -1);
    1863        5262 :   GEN a = ncV_polint_center_tree(A, P, T, R, m2);
    1864        5262 :   return gerepileupto(av, a);
    1865             : }
    1866             : 
    1867             : GEN
    1868           0 : nmV_chinese_center_tree(GEN A, GEN P, GEN T, GEN R)
    1869             : {
    1870           0 :   pari_sp av = avma;
    1871           0 :   GEN m2 = shifti(gmael(T, lg(T)-1, 1), -1);
    1872           0 :   GEN a = nmV_polint_center_tree(A, P, T, R, m2);
    1873           0 :   return gerepileupto(av, a);
    1874             : }
    1875             : 
    1876             : GEN
    1877       52975 : nmV_chinese_center_tree_seq(GEN A, GEN P, GEN T, GEN R)
    1878             : {
    1879       52975 :   pari_sp av = avma;
    1880       52975 :   GEN m2 = shifti(gmael(T, lg(T)-1, 1), -1);
    1881       52973 :   GEN a = nmV_polint_center_tree_seq(A, P, T, R, m2);
    1882       52976 :   return gerepileupto(av, a);
    1883             : }
    1884             : 
    1885             : GEN
    1886      114543 : nmV_chinese_center(GEN A, GEN P, GEN *pt_mod)
    1887             : {
    1888      114543 :   pari_sp av = avma;
    1889      114543 :   GEN T = ZV_producttree(P);
    1890      114543 :   GEN R = ZV_chinesetree(P, T);
    1891      114543 :   GEN m2 = shifti(gmael(T, lg(T)-1, 1), -1);
    1892      114543 :   GEN a = nmV_polint_center_tree(A, P, T, R, m2);
    1893      114543 :   return gc_chinese(av, T, a, pt_mod);
    1894             : }
    1895             : 
    1896             : GEN
    1897           0 : nxCV_chinese_center_tree(GEN A, GEN P, GEN T, GEN R)
    1898             : {
    1899           0 :   pari_sp av = avma;
    1900           0 :   GEN m2 = shifti(gmael(T, lg(T)-1, 1), -1);
    1901           0 :   GEN a = nxCV_polint_center_tree(A, P, T, R, m2);
    1902           0 :   return gerepileupto(av, a);
    1903             : }
    1904             : 
    1905             : GEN
    1906           0 : nxCV_chinese_center(GEN A, GEN P, GEN *pt_mod)
    1907             : {
    1908           0 :   pari_sp av = avma;
    1909           0 :   GEN T = ZV_producttree(P);
    1910           0 :   GEN R = ZV_chinesetree(P, T);
    1911           0 :   GEN m2 = shifti(gmael(T, lg(T)-1, 1), -1);
    1912           0 :   GEN a = nxCV_polint_center_tree(A, P, T, R, m2);
    1913           0 :   return gc_chinese(av, T, a, pt_mod);
    1914             : }
    1915             : 
    1916             : GEN
    1917         430 : nxMV_chinese_center_tree_seq(GEN A, GEN P, GEN T, GEN R)
    1918             : {
    1919         430 :   pari_sp av = avma;
    1920         430 :   GEN m2 = shifti(gmael(T, lg(T)-1, 1), -1);
    1921         430 :   GEN a = nxMV_polint_center_tree_seq(A, P, T, R, m2);
    1922         430 :   return gerepileupto(av, a);
    1923             : }
    1924             : 
    1925             : GEN
    1926          90 : nxMV_chinese_center(GEN A, GEN P, GEN *pt_mod)
    1927             : {
    1928          90 :   pari_sp av = avma;
    1929          90 :   GEN T = ZV_producttree(P);
    1930          90 :   GEN R = ZV_chinesetree(P, T);
    1931          90 :   GEN m2 = shifti(gmael(T, lg(T)-1, 1), -1);
    1932          90 :   GEN a = nxMV_polint_center_tree(A, P, T, R, m2);
    1933          90 :   return gc_chinese(av, T, a, pt_mod);
    1934             : }
    1935             : 
    1936             : /**********************************************************************
    1937             :  **                    Powering  over (Z/NZ)^*, small N              **
    1938             :  **********************************************************************/
    1939             : 
    1940             : /* 2^n mod p; assume n > 1 */
    1941             : static ulong
    1942    11576477 : Fl_2powu_pre(ulong n, ulong p, ulong pi)
    1943             : {
    1944    11576477 :   ulong y = 2;
    1945    11576477 :   int j = 1+bfffo(n);
    1946             :   /* normalize, i.e set highest bit to 1 (we know n != 0) */
    1947    11576477 :   n<<=j; j = BITS_IN_LONG-j; /* first bit is now implicit */
    1948   390056122 :   for (; j; n<<=1,j--)
    1949             :   {
    1950   378480014 :     y = Fl_sqr_pre(y,p,pi);
    1951   378449766 :     if (n & HIGHBIT) y = Fl_double(y, p);
    1952             :   }
    1953    11576108 :   return y;
    1954             : }
    1955             : 
    1956             : /* 2^n mod p; assume n > 1 and !(p & HIGHMASK) */
    1957             : static ulong
    1958     1452803 : Fl_2powu(ulong n, ulong p)
    1959             : {
    1960     1452803 :   ulong y = 2;
    1961     1452803 :   int j = 1+bfffo(n);
    1962             :   /* normalize, i.e set highest bit to 1 (we know n != 0) */
    1963     1452803 :   n<<=j; j = BITS_IN_LONG-j; /* first bit is now implicit */
    1964     5597052 :   for (; j; n<<=1,j--)
    1965             :   {
    1966     4144240 :     y = (y*y) % p;
    1967     4144240 :     if (n & HIGHBIT) y = Fl_double(y, p);
    1968             :   }
    1969     1452812 :   return y;
    1970             : }
    1971             : 
    1972             : ulong
    1973    52218133 : Fl_powu_pre(ulong x, ulong n0, ulong p, ulong pi)
    1974             : {
    1975             :   ulong y, z, n;
    1976    52218133 :   if (n0 <= 1)
    1977             :   { /* frequent special cases */
    1978     5063779 :     if (n0 == 1) return x;
    1979     1346211 :     if (n0 == 0) return 1;
    1980             :   }
    1981    47154293 :   if (x <= 2)
    1982             :   {
    1983    12124379 :     if (x == 2) return Fl_2powu_pre(n0, p, pi);
    1984      546823 :     return x; /* 0 or 1 */
    1985             :   }
    1986    35029914 :   y = 1; z = x; n = n0;
    1987             :   for(;;)
    1988             :   {
    1989   448381005 :     if (n&1) y = Fl_mul_pre(y,z,p,pi);
    1990   448556944 :     n>>=1; if (!n) return y;
    1991   413556304 :     z = Fl_sqr_pre(z,p,pi);
    1992             :   }
    1993             : }
    1994             : 
    1995             : ulong
    1996   147114994 : Fl_powu(ulong x, ulong n0, ulong p)
    1997             : {
    1998             :   ulong y, z, n;
    1999   147114994 :   if (n0 <= 2)
    2000             :   { /* frequent special cases */
    2001    92252795 :     if (n0 == 2) return Fl_sqr(x,p);
    2002    23395317 :     if (n0 == 1) return x;
    2003      500151 :     if (n0 == 0) return 1;
    2004             :   }
    2005    54820297 :   if (x <= 1) return x; /* 0 or 1 */
    2006    54605532 :   if (p & HIGHMASK)
    2007     6317281 :     return Fl_powu_pre(x, n0, p, get_Fl_red(p));
    2008    48288251 :   if (x == 2) return Fl_2powu(n0, p);
    2009    46835447 :   y = 1; z = x; n = n0;
    2010             :   for(;;)
    2011             :   {
    2012   253193491 :     if (n&1) y = (y*z) % p;
    2013   253193491 :     n>>=1; if (!n) return y;
    2014   206358044 :     z = (z*z) % p;
    2015             :   }
    2016             : }
    2017             : 
    2018             : /* Reduce data dependency to maximize internal parallelism */
    2019             : GEN
    2020    11330899 : Fl_powers_pre(ulong x, long n, ulong p, ulong pi)
    2021             : {
    2022             :   long i, k;
    2023    11330899 :   GEN powers = cgetg(n + 2, t_VECSMALL);
    2024    11318671 :   powers[1] = 1; if (n == 0) return powers;
    2025    11318671 :   powers[2] = x;
    2026   282845272 :   for (i = 3, k=2; i <= n; i+=2, k++)
    2027             :   {
    2028   271510444 :     powers[i] = Fl_sqr_pre(powers[k], p, pi);
    2029   271511209 :     powers[i+1] = Fl_mul_pre(powers[k], powers[k+1], p, pi);
    2030             :   }
    2031    11334828 :   if (i==n+1)
    2032     9799232 :     powers[i] = Fl_sqr_pre(powers[k], p, pi);
    2033    11335632 :   return powers;
    2034             : }
    2035             : 
    2036             : GEN
    2037       24981 : Fl_powers(ulong x, long n, ulong p)
    2038             : {
    2039       24981 :   return Fl_powers_pre(x, n, p, get_Fl_red(p));
    2040             : }
    2041             : 
    2042             : /**********************************************************************
    2043             :  **                    Powering  over (Z/NZ)^*, large N              **
    2044             :  **********************************************************************/
    2045             : 
    2046             : static GEN
    2047     4708642 : Fp_dblsqr(GEN x, GEN N)
    2048             : {
    2049     4708642 :   GEN z = shifti(Fp_sqr(x, N), 1);
    2050     4587246 :   return cmpii(z, N) >= 0? subii(z, N): z;
    2051             : }
    2052             : 
    2053             : typedef struct muldata {
    2054             :   GEN (*sqr)(void * E, GEN x);
    2055             :   GEN (*mul)(void * E, GEN x, GEN y);
    2056             :   GEN (*mul2)(void * E, GEN x);
    2057             : } muldata;
    2058             : 
    2059             : /* modified Barrett reduction with one fold */
    2060             : /* See Fast Modular Reduction, W. Hasenplaugh, G. Gaubatz, V. Gopal, ARITH 18 */
    2061             : 
    2062             : static GEN
    2063        8611 : Fp_invmBarrett(GEN p, long s)
    2064             : {
    2065        8611 :   GEN R, Q = dvmdii(int2n(3*s),p,&R);
    2066        8611 :   return mkvec2(Q,R);
    2067             : }
    2068             : 
    2069             : /* a <= (N-1)^2, 2^(2s-2) <= N < 2^(2s). Return 0 <= r < N such that
    2070             :  * a = r (mod N) */
    2071             : static GEN
    2072     4698940 : Fp_rem_mBarrett(GEN a, GEN B, long s, GEN N)
    2073             : {
    2074     4698940 :   pari_sp av = avma;
    2075     4698940 :   GEN P = gel(B, 1), Q = gel(B, 2); /* 2^(3s) = P N + Q, 0 <= Q < N */
    2076     4698940 :   long t = expi(P)+1; /* 2^(t-1) <= P < 2^t */
    2077     4698940 :   GEN u = shifti(a, -3*s), v = remi2n(a, 3*s); /* a = 2^(3s)u + v */
    2078     4698940 :   GEN A = addii(v, mulii(Q,u)); /* 0 <= A < 2^(3s+1) */
    2079     4698940 :   GEN q = shifti(mulii(shifti(A, t-3*s), P), -t); /* A/N - 4 < q <= A/N */
    2080     4698940 :   GEN r = subii(A, mulii(q, N));
    2081     4698940 :   GEN sr= subii(r,N);     /* 0 <= r < 4*N */
    2082     4698940 :   if (signe(sr)<0) return gerepileuptoint(av, r);
    2083     2787452 :   r=sr; sr = subii(r,N);  /* 0 <= r < 3*N */
    2084     2787452 :   if (signe(sr)<0) return gerepileuptoint(av, r);
    2085      107472 :   r=sr; sr = subii(r,N);  /* 0 <= r < 2*N */
    2086      107472 :   return gerepileuptoint(av, signe(sr)>=0 ? sr:r);
    2087             : }
    2088             : 
    2089             : /* Montgomery reduction */
    2090             : 
    2091             : INLINE ulong
    2092      807181 : init_montdata(GEN N) { return (ulong) -invmod2BIL(mod2BIL(N)); }
    2093             : 
    2094             : struct montred
    2095             : {
    2096             :   GEN N;
    2097             :   ulong inv;
    2098             : };
    2099             : 
    2100             : /* Montgomery reduction */
    2101             : static GEN
    2102    49771465 : _sqr_montred(void * E, GEN x)
    2103             : {
    2104    49771465 :   struct montred * D = (struct montred *) E;
    2105    49771465 :   return red_montgomery(sqri(x), D->N, D->inv);
    2106             : }
    2107             : 
    2108             : /* Montgomery reduction */
    2109             : static GEN
    2110     4656997 : _mul_montred(void * E, GEN x, GEN y)
    2111             : {
    2112     4656997 :   struct montred * D = (struct montred *) E;
    2113     4656997 :   return red_montgomery(mulii(x, y), D->N, D->inv);
    2114             : }
    2115             : 
    2116             : static GEN
    2117     7528967 : _mul2_montred(void * E, GEN x)
    2118             : {
    2119     7528967 :   struct montred * D = (struct montred *) E;
    2120     7528967 :   GEN z = shifti(_sqr_montred(E, x), 1);
    2121     7526257 :   long l = lgefint(D->N);
    2122     7969358 :   while (lgefint(z) > l) z = subii(z, D->N);
    2123     7526832 :   return z;
    2124             : }
    2125             : 
    2126             : static GEN
    2127    19737303 : _sqr_remii(void* N, GEN x)
    2128    19737303 : { return remii(sqri(x), (GEN) N); }
    2129             : 
    2130             : static GEN
    2131     1432600 : _mul_remii(void* N, GEN x, GEN y)
    2132     1432600 : { return remii(mulii(x, y), (GEN) N); }
    2133             : 
    2134             : static GEN
    2135     3292437 : _mul2_remii(void* N, GEN x)
    2136     3292437 : { return Fp_dblsqr(x, (GEN) N); }
    2137             : 
    2138             : struct redbarrett
    2139             : {
    2140             :   GEN iM, N;
    2141             :   long s;
    2142             : };
    2143             : 
    2144             : static GEN
    2145     4243321 : _sqr_remiibar(void *E, GEN x)
    2146             : {
    2147     4243321 :   struct redbarrett * D = (struct redbarrett *) E;
    2148     4243321 :   return Fp_rem_mBarrett(sqri(x), D->iM, D->s, D->N);
    2149             : }
    2150             : 
    2151             : static GEN
    2152      455619 : _mul_remiibar(void *E, GEN x, GEN y)
    2153             : {
    2154      455619 :   struct redbarrett * D = (struct redbarrett *) E;
    2155      455619 :   return Fp_rem_mBarrett(mulii(x, y), D->iM, D->s, D->N);
    2156             : }
    2157             : 
    2158             : static GEN
    2159     1413630 : _mul2_remiibar(void *E, GEN x)
    2160             : {
    2161     1413630 :   struct redbarrett * D = (struct redbarrett *) E;
    2162     1413630 :   return Fp_dblsqr(x, D->N);
    2163             : }
    2164             : 
    2165             : static long
    2166     1017199 : Fp_select_red(GEN *y, ulong k, GEN N, long lN, muldata *D, void **pt_E)
    2167             : {
    2168     1017199 :   if (lN >= Fp_POW_BARRETT_LIMIT && (k==0 || ((double)k)*expi(*y) > 2 + expi(N)))
    2169             :   {
    2170        8611 :     struct redbarrett * E = (struct redbarrett *) stack_malloc(sizeof(struct redbarrett));
    2171        8611 :     D->sqr = &_sqr_remiibar;
    2172        8611 :     D->mul = &_mul_remiibar;
    2173        8611 :     D->mul2 = &_mul2_remiibar;
    2174        8611 :     E->N = N;
    2175        8611 :     E->s = 1+(expi(N)>>1);
    2176        8611 :     E->iM = Fp_invmBarrett(N, E->s);
    2177        8611 :     *pt_E = (void*) E;
    2178        8611 :     return 0;
    2179             :   }
    2180     1008588 :   else if (mod2(N) && lN < Fp_POW_REDC_LIMIT)
    2181             :   {
    2182      807190 :     struct montred * E = (struct montred *) stack_malloc(sizeof(struct montred));
    2183      807190 :     *y = remii(shifti(*y, bit_accuracy(lN)), N);
    2184      807185 :     D->sqr = &_sqr_montred;
    2185      807185 :     D->mul = &_mul_montred;
    2186      807185 :     D->mul2 = &_mul2_montred;
    2187      807185 :     E->N = N;
    2188      807185 :     E->inv = init_montdata(N);
    2189      807175 :     *pt_E = (void*) E;
    2190      807175 :     return 1;
    2191             :   }
    2192             :   else
    2193             :   {
    2194      201393 :     D->sqr = &_sqr_remii;
    2195      201393 :     D->mul = &_mul_remii;
    2196      201393 :     D->mul2 = &_mul2_remii;
    2197      201393 :     *pt_E = (void*) N;
    2198      201393 :     return 0;
    2199             :   }
    2200             : }
    2201             : 
    2202             : GEN
    2203     3055308 : Fp_powu(GEN A, ulong k, GEN N)
    2204             : {
    2205     3055308 :   long lN = lgefint(N);
    2206             :   int base_is_2, use_montgomery;
    2207             :   muldata D;
    2208             :   void *E;
    2209             :   pari_sp av;
    2210             : 
    2211     3055308 :   if (lN == 3) {
    2212     1432839 :     ulong n = uel(N,2);
    2213     1432839 :     return utoi( Fl_powu(umodiu(A, n), k, n) );
    2214             :   }
    2215     1622469 :   if (k <= 2)
    2216             :   { /* frequent special cases */
    2217      814860 :     if (k == 2) return Fp_sqr(A,N);
    2218      275650 :     if (k == 1) return A;
    2219           0 :     if (k == 0) return gen_1;
    2220             :   }
    2221      807608 :   av = avma; A = modii(A,N);
    2222      807609 :   base_is_2 = 0;
    2223      807609 :   if (lgefint(A) == 3) switch(A[2])
    2224             :   {
    2225         770 :     case 1: set_avma(av); return gen_1;
    2226       38967 :     case 2:  base_is_2 = 1; break;
    2227             :   }
    2228             : 
    2229             :   /* TODO: Move this out of here and use for general modular computations */
    2230      806839 :   use_montgomery = Fp_select_red(&A, k, N, lN, &D, &E);
    2231      806839 :   if (base_is_2)
    2232       38967 :     A = gen_powu_fold_i(A, k, E, D.sqr, D.mul2);
    2233             :   else
    2234      767872 :     A = gen_powu_i(A, k, E, D.sqr, D.mul);
    2235      806839 :   if (use_montgomery)
    2236             :   {
    2237      686756 :     A = red_montgomery(A, N, ((struct montred *) E)->inv);
    2238      686756 :     if (cmpii(A, N) >= 0) A = subii(A,N);
    2239             :   }
    2240      806839 :   return gerepileuptoint(av, A);
    2241             : }
    2242             : 
    2243             : GEN
    2244       22428 : Fp_pows(GEN A, long k, GEN N)
    2245             : {
    2246       22428 :   if (lgefint(N) == 3) {
    2247        7939 :     ulong n = N[2];
    2248        7939 :     ulong a = umodiu(A, n);
    2249        7939 :     if (k < 0) {
    2250         133 :       a = Fl_inv(a, n);
    2251         133 :       k = -k;
    2252             :     }
    2253        7939 :     return utoi( Fl_powu(a, (ulong)k, n) );
    2254             :   }
    2255       14489 :   if (k < 0) { A = Fp_inv(A, N); k = -k; };
    2256       14489 :   return Fp_powu(A, (ulong)k, N);
    2257             : }
    2258             : 
    2259             : /* A^K mod N */
    2260             : GEN
    2261    34430762 : Fp_pow(GEN A, GEN K, GEN N)
    2262             : {
    2263             :   pari_sp av;
    2264    34430762 :   long s, lN = lgefint(N), sA, sy;
    2265             :   int base_is_2, use_montgomery;
    2266             :   GEN y;
    2267             :   muldata D;
    2268             :   void *E;
    2269             : 
    2270    34430762 :   s = signe(K);
    2271    34430762 :   if (!s) return dvdii(A,N)? gen_0: gen_1;
    2272    33401821 :   if (lN == 3 && lgefint(K) == 3)
    2273             :   {
    2274    32508094 :     ulong n = N[2], a = umodiu(A, n);
    2275    32508451 :     if (s < 0) a = Fl_inv(a, n);
    2276    32508520 :     if (a <= 1) return utoi(a); /* 0 or 1 */
    2277    29547888 :     return utoi(Fl_powu(a, uel(K,2), n));
    2278             :   }
    2279             : 
    2280      893727 :   av = avma;
    2281      893727 :   if (s < 0) y = Fp_inv(A,N);
    2282             :   else
    2283             :   {
    2284      891770 :     y = modii(A,N);
    2285      892008 :     if (!signe(y)) { set_avma(av); return gen_0; }
    2286             :   }
    2287      893966 :   if (lgefint(K) == 3) return gerepileuptoint(av, Fp_powu(y, K[2], N));
    2288             : 
    2289      210536 :   base_is_2 = 0;
    2290      210536 :   sy = abscmpii(y, shifti(N,-1)) > 0;
    2291      210561 :   if (sy) y = subii(N,y);
    2292      210565 :   sA = sy && mod2(K);
    2293      210563 :   if (lgefint(y) == 3) switch(y[2])
    2294             :   {
    2295         199 :     case 1:  set_avma(av); return sA ? subis(N,1): gen_1;
    2296      139823 :     case 2:  base_is_2 = 1; break;
    2297             :   }
    2298             : 
    2299             :   /* TODO: Move this out of here and use for general modular computations */
    2300      210364 :   use_montgomery = Fp_select_red(&y, 0UL, N, lN, &D, &E);
    2301      210348 :   if (base_is_2)
    2302      139833 :     y = gen_pow_fold_i(y, K, E, D.sqr, D.mul2);
    2303             :   else
    2304       70515 :     y = gen_pow_i(y, K, E, D.sqr, D.mul);
    2305      210376 :   if (use_montgomery)
    2306             :   {
    2307      120442 :     y = red_montgomery(y, N, ((struct montred *) E)->inv);
    2308      120443 :     if (cmpii(y,N) >= 0) y = subii(y,N);
    2309             :   }
    2310      210376 :   if (sA) y = subii(N, y);
    2311      210376 :   return gerepileuptoint(av,y);
    2312             : }
    2313             : 
    2314             : static GEN
    2315    14030630 : _Fp_mul(void *E, GEN x, GEN y) { return Fp_mul(x,y,(GEN)E); }
    2316             : static GEN
    2317     8126876 : _Fp_sqr(void *E, GEN x) { return Fp_sqr(x,(GEN)E); }
    2318             : static GEN
    2319       56163 : _Fp_one(void *E) { (void) E; return gen_1; }
    2320             : 
    2321             : GEN
    2322          77 : Fp_pow_init(GEN x, GEN n, long k, GEN p)
    2323          77 : { return gen_pow_init(x, n, k, (void*)p, &_Fp_sqr, &_Fp_mul); }
    2324             : 
    2325             : GEN
    2326       54096 : Fp_pow_table(GEN R, GEN n, GEN p)
    2327       54096 : { return gen_pow_table(R, n, (void*)p, &_Fp_one, &_Fp_mul); }
    2328             : 
    2329             : GEN
    2330        4523 : Fp_powers(GEN x, long n, GEN p)
    2331             : {
    2332        4523 :   if (lgefint(p) == 3)
    2333        2456 :     return Flv_to_ZV(Fl_powers(umodiu(x, uel(p, 2)), n, uel(p, 2)));
    2334        2067 :   return gen_powers(x, n, 1, (void*)p, _Fp_sqr, _Fp_mul, _Fp_one);
    2335             : }
    2336             : 
    2337             : GEN
    2338         434 : FpV_prod(GEN V, GEN p) { return gen_product(V, (void *)p, &_Fp_mul); }
    2339             : 
    2340             : static GEN
    2341    27110731 : _Fp_pow(void *E, GEN x, GEN n) { return Fp_pow(x,n,(GEN)E); }
    2342             : static GEN
    2343         147 : _Fp_rand(void *E) { return addiu(randomi(subiu((GEN)E,1)),1); }
    2344             : 
    2345             : static GEN Fp_easylog(void *E, GEN a, GEN g, GEN ord);
    2346             : static const struct bb_group Fp_star={_Fp_mul,_Fp_pow,_Fp_rand,hash_GEN,
    2347             :                                       equalii,equali1,Fp_easylog};
    2348             : 
    2349             : static GEN
    2350      761553 : _Fp_red(void *E, GEN x) { return Fp_red(x, (GEN)E); }
    2351             : static GEN
    2352      853646 : _Fp_add(void *E, GEN x, GEN y) { (void) E; return addii(x,y); }
    2353             : static GEN
    2354      767753 : _Fp_neg(void *E, GEN x) { (void) E; return negi(x); }
    2355             : static GEN
    2356      497784 : _Fp_rmul(void *E, GEN x, GEN y) { (void) E; return mulii(x,y); }
    2357             : static GEN
    2358       33090 : _Fp_inv(void *E, GEN x) { return Fp_inv(x,(GEN)E); }
    2359             : static int
    2360      218619 : _Fp_equal0(GEN x) { return signe(x)==0; }
    2361             : static GEN
    2362       28761 : _Fp_s(void *E, long x) { (void) E; return stoi(x); }
    2363             : 
    2364             : static const struct bb_field Fp_field={_Fp_red,_Fp_add,_Fp_rmul,_Fp_neg,
    2365             :                                         _Fp_inv,_Fp_equal0,_Fp_s};
    2366             : 
    2367        7077 : const struct bb_field *get_Fp_field(void **E, GEN p)
    2368        7077 : { *E = (void*)p; return &Fp_field; }
    2369             : 
    2370             : /*********************************************************************/
    2371             : /**               ORDER of INTEGERMOD x  in  (Z/nZ)*                **/
    2372             : /*********************************************************************/
    2373             : ulong
    2374       15078 : Fl_order(ulong a, ulong o, ulong p)
    2375             : {
    2376       15078 :   pari_sp av = avma;
    2377             :   GEN m, P, E;
    2378             :   long i;
    2379       15078 :   if (a==1) return 1;
    2380        9954 :   if (!o) o = p-1;
    2381        9954 :   m = factoru(o);
    2382        9954 :   P = gel(m,1);
    2383        9954 :   E = gel(m,2);
    2384       25018 :   for (i = lg(P)-1; i; i--)
    2385             :   {
    2386       15064 :     ulong j, l = P[i], e = E[i], t = o / upowuu(l,e), y = Fl_powu(a, t, p);
    2387       15064 :     if (y == 1) o = t;
    2388       16779 :     else for (j = 1; j < e; j++)
    2389             :     {
    2390        5257 :       y = Fl_powu(y, l, p);
    2391        5257 :       if (y == 1) { o = t *  upowuu(l, j); break; }
    2392             :     }
    2393             :   }
    2394        9954 :   return gc_ulong(av, o);
    2395             : }
    2396             : 
    2397             : /*Find the exact order of a assuming a^o==1*/
    2398             : GEN
    2399       10438 : Fp_order(GEN a, GEN o, GEN p) {
    2400       10438 :   if (lgefint(p) == 3 && (!o || typ(o) == t_INT))
    2401             :   {
    2402          21 :     ulong pp = p[2], oo = (o && lgefint(o)==3)? uel(o,2): pp-1;
    2403          21 :     return utoi( Fl_order(umodiu(a, pp), oo, pp) );
    2404             :   }
    2405       10417 :   return gen_order(a, o, (void*)p, &Fp_star);
    2406             : }
    2407             : GEN
    2408          70 : Fp_factored_order(GEN a, GEN o, GEN p)
    2409          70 : { return gen_factored_order(a, o, (void*)p, &Fp_star); }
    2410             : 
    2411             : /* return order of a mod p^e, e > 0, pe = p^e */
    2412             : static GEN
    2413          70 : Zp_order(GEN a, GEN p, long e, GEN pe)
    2414             : {
    2415             :   GEN ap, op;
    2416          70 :   if (absequaliu(p, 2))
    2417             :   {
    2418          56 :     if (e == 1) return gen_1;
    2419          56 :     if (e == 2) return mod4(a) == 1? gen_1: gen_2;
    2420          49 :     if (mod4(a) == 1) op = gen_1; else { op = gen_2; a = Fp_sqr(a, pe); }
    2421             :   } else {
    2422          14 :     ap = (e == 1)? a: remii(a,p);
    2423          14 :     op = Fp_order(ap, subiu(p,1), p);
    2424          14 :     if (e == 1) return op;
    2425           0 :     a = Fp_pow(a, op, pe); /* 1 mod p */
    2426             :   }
    2427          49 :   if (equali1(a)) return op;
    2428           7 :   return mulii(op, powiu(p, e - Z_pval(subiu(a,1), p)));
    2429             : }
    2430             : 
    2431             : GEN
    2432          63 : znorder(GEN x, GEN o)
    2433             : {
    2434          63 :   pari_sp av = avma;
    2435             :   GEN b, a;
    2436             : 
    2437          63 :   if (typ(x) != t_INTMOD) pari_err_TYPE("znorder [t_INTMOD expected]",x);
    2438          56 :   b = gel(x,1); a = gel(x,2);
    2439          56 :   if (!equali1(gcdii(a,b))) pari_err_COPRIME("znorder", a,b);
    2440          49 :   if (!o)
    2441             :   {
    2442          35 :     GEN fa = Z_factor(b), P = gel(fa,1), E = gel(fa,2);
    2443          35 :     long i, l = lg(P);
    2444          35 :     o = gen_1;
    2445          70 :     for (i = 1; i < l; i++)
    2446             :     {
    2447          35 :       GEN p = gel(P,i);
    2448          35 :       long e = itos(gel(E,i));
    2449             : 
    2450          35 :       if (l == 2)
    2451          35 :         o = Zp_order(a, p, e, b);
    2452             :       else {
    2453           0 :         GEN pe = powiu(p,e);
    2454           0 :         o = lcmii(o, Zp_order(remii(a,pe), p, e, pe));
    2455             :       }
    2456             :     }
    2457          35 :     return gerepileuptoint(av, o);
    2458             :   }
    2459          14 :   return Fp_order(a, o, b);
    2460             : }
    2461             : 
    2462             : /*********************************************************************/
    2463             : /**               DISCRETE LOGARITHM  in  (Z/nZ)*                   **/
    2464             : /*********************************************************************/
    2465             : static GEN
    2466       58614 : Fp_log_halfgcd(ulong bnd, GEN C, GEN g, GEN p)
    2467             : {
    2468       58614 :   pari_sp av = avma;
    2469             :   GEN h1, h2, F, G;
    2470       58614 :   if (!Fp_ratlift(g,p,C,shifti(C,-1),&h1,&h2)) return gc_NULL(av);
    2471       35205 :   if ((F = Z_issmooth_fact(h1, bnd)) && (G = Z_issmooth_fact(h2, bnd)))
    2472             :   {
    2473         126 :     GEN M = cgetg(3, t_MAT);
    2474         126 :     gel(M,1) = vecsmall_concat(gel(F, 1),gel(G, 1));
    2475         126 :     gel(M,2) = vecsmall_concat(gel(F, 2),zv_neg_inplace(gel(G, 2)));
    2476         126 :     return gerepileupto(av, M);
    2477             :   }
    2478       35079 :   return gc_NULL(av);
    2479             : }
    2480             : 
    2481             : static GEN
    2482       58614 : Fp_log_find_rel(GEN b, ulong bnd, GEN C, GEN p, GEN *g, long *e)
    2483             : {
    2484             :   GEN rel;
    2485       58614 :   do { (*e)++; *g = Fp_mul(*g, b, p); rel = Fp_log_halfgcd(bnd, C, *g, p); }
    2486       58614 :   while (!rel);
    2487         126 :   return rel;
    2488             : }
    2489             : 
    2490             : struct Fp_log_rel
    2491             : {
    2492             :   GEN rel;
    2493             :   ulong prmax;
    2494             :   long nbrel, nbmax, nbgen;
    2495             : };
    2496             : 
    2497             : /* add u^e */
    2498             : static void
    2499        3157 : addifsmooth1(struct Fp_log_rel *r, GEN z, long u, long e)
    2500             : {
    2501        3157 :   pari_sp av = avma;
    2502        3157 :   long off = r->prmax+1;
    2503        3157 :   GEN F = cgetg(3, t_MAT);
    2504        3157 :   gel(F,1) = vecsmall_append(gel(z,1), off+u);
    2505        3157 :   gel(F,2) = vecsmall_append(gel(z,2), e);
    2506        3157 :   gel(r->rel,++r->nbrel) = gerepileupto(av, F);
    2507        3157 : }
    2508             : 
    2509             : /* add u^-1 v^-1 */
    2510             : static void
    2511      104083 : addifsmooth2(struct Fp_log_rel *r, GEN z, long u, long v)
    2512             : {
    2513      104083 :   pari_sp av = avma;
    2514      104083 :   long off = r->prmax+1;
    2515      104083 :   GEN P = mkvecsmall2(off+u,off+v), E = mkvecsmall2(-1,-1);
    2516      104083 :   GEN F = cgetg(3, t_MAT);
    2517      104083 :   gel(F,1) = vecsmall_concat(gel(z,1), P);
    2518      104083 :   gel(F,2) = vecsmall_concat(gel(z,2), E);
    2519      104083 :   gel(r->rel,++r->nbrel) = gerepileupto(av, F);
    2520      104083 : }
    2521             : 
    2522             : /*
    2523             : Let p=C^2+c
    2524             : Solve h = (C+x)*(C+a)-p = 0 [mod l]
    2525             : h= -c+x*(C+a)+C*a = 0  [mod l]
    2526             : x = (c-C*a)/(C+a) [mod l]
    2527             : h = -c+C*(x+a)+a*x
    2528             : */
    2529             : 
    2530             : GEN
    2531       40827 : Fp_log_sieve_worker(long a, long prmax, GEN C, GEN c, GEN Ci, GEN ci, GEN pi, GEN sz)
    2532             : {
    2533       40827 :   pari_sp ltop = avma;
    2534       40827 :   long i, j, th, n = lg(pi)-1, rel = 1;
    2535       40827 :   GEN sieve = zero_zv(a+2)+1;
    2536       40850 :   GEN L = cgetg(1+a+2, t_VEC);
    2537       40837 :   pari_sp av = avma;
    2538       40837 :   GEN z, h = addis(C,a);
    2539       40817 :   if ((z = Z_issmooth_fact(h, prmax)))
    2540             :   {
    2541        3007 :     gel(L, rel++) = mkvec2(z, mkvecsmall3(1, a, -1));
    2542        3016 :     av = avma;
    2543             :   }
    2544    16835364 :   for (i=1; i<=n; i++)
    2545             :   {
    2546    16809850 :     ulong li = pi[i], s = sz[i], al = a % li;
    2547    16809850 :     ulong u, iv = Fl_invsafe(Fl_add(Ci[i],al,li),li);
    2548    17154451 :     if (!iv) continue;
    2549    16730061 :     u = Fl_mul(Fl_sub(ci[i],Fl_mul(Ci[i],al,li),li), iv ,li);
    2550    77389417 :     for(j = u; j<=a; j+=li) sieve[j] += s;
    2551             :   }
    2552       35046 :   if (a)
    2553             :   {
    2554       40737 :     long e = expi(mulis(C,a));
    2555       40798 :     th = e - expu(e) - 1;
    2556          54 :   } else th = -1;
    2557    28023325 :   for (j=0; j<a; j++)
    2558    27981761 :     if (sieve[j]>=th)
    2559             :     {
    2560      119453 :       GEN h = addiu(subii(muliu(C,a+j),c), a*j);
    2561      119334 :       if ((z = Z_issmooth_fact(h, prmax)))
    2562             :       {
    2563      109798 :         gel(L, rel++) = mkvec2(z, mkvecsmall3(2, a, j));
    2564      109993 :         av = avma;
    2565        9431 :       } else set_avma(av);
    2566             :     }
    2567             :   /* j = a */
    2568       41564 :   if (sieve[a]>=th)
    2569             :   {
    2570         476 :     GEN h = addiu(subii(muliu(C,2*a),c), a*a);
    2571         476 :     if ((z = Z_issmooth_fact(h, prmax)))
    2572         385 :       gel(L, rel++) = mkvec2(z, mkvecsmall3(1, a, -2));
    2573             :   }
    2574       41564 :   setlg(L, rel); return gerepilecopy(ltop, L);
    2575             : }
    2576             : 
    2577             : static long
    2578          63 : Fp_log_sieve(struct Fp_log_rel *r, GEN C, GEN c, GEN Ci, GEN ci, GEN pi, GEN sz)
    2579             : {
    2580             :   struct pari_mt pt;
    2581          63 :   long i, j, nb = 0;
    2582          63 :   GEN worker = snm_closure(is_entry("_Fp_log_sieve_worker"),
    2583             :                mkvecn(7, utoi(r->prmax), C, c, Ci, ci, pi, sz));
    2584          63 :   long running, pending = 0;
    2585          63 :   GEN W = zerovec(r->nbgen);
    2586          63 :   mt_queue_start_lim(&pt, worker, r->nbgen);
    2587       41229 :   for (i = 0; (running = (i < r->nbgen)) || pending; i++)
    2588             :   {
    2589             :     GEN done;
    2590             :     long idx;
    2591       41166 :     mt_queue_submit(&pt, i, running ? mkvec(stoi(i)): NULL);
    2592       41166 :     done = mt_queue_get(&pt, &idx, &pending);
    2593       41166 :     if (!done || lg(done)==1) continue;
    2594       35917 :     gel(W, idx+1) = done;
    2595       35917 :     nb += lg(done)-1;
    2596       35917 :     if (DEBUGLEVEL && (i&127)==0)
    2597           0 :       err_printf("%ld%% ",100*nb/r->nbmax);
    2598             :   }
    2599          63 :   mt_queue_end(&pt);
    2600       39550 :   for(j = 1; j <= r->nbgen && r->nbrel < r->nbmax; j++)
    2601             :   {
    2602             :     long ll, m;
    2603       39487 :     GEN L = gel(W,j);
    2604       39487 :     if (isintzero(L)) continue;
    2605       34531 :     ll = lg(L);
    2606      141771 :     for (m=1; m<ll && r->nbrel < r->nbmax ; m++)
    2607             :     {
    2608      107240 :       GEN Lm = gel(L,m), h = gel(Lm, 1), v = gel(Lm, 2);
    2609      107240 :       if (v[1] == 1)
    2610        3157 :         addifsmooth1(r, h, v[2], v[3]);
    2611             :       else
    2612      104083 :         addifsmooth2(r, h, v[2], v[3]);
    2613             :     }
    2614             :   }
    2615          63 :   return j;
    2616             : }
    2617             : 
    2618             : static GEN
    2619         735 : ECP_psi(GEN x, GEN y)
    2620             : {
    2621         735 :   long prec = realprec(x);
    2622         735 :   GEN lx = glog(x, prec), ly = glog(y, prec);
    2623         735 :   GEN u = gdiv(lx, ly);
    2624         735 :   return gpow(u, gneg(u),prec);
    2625             : }
    2626             : 
    2627             : struct computeG
    2628             : {
    2629             :   GEN C;
    2630             :   long bnd, nbi;
    2631             : };
    2632             : 
    2633             : static GEN
    2634         735 : _computeG(void *E, GEN gen)
    2635             : {
    2636         735 :   struct computeG * d = (struct computeG *) E;
    2637         735 :   GEN ps = ECP_psi(gmul(gen,d->C), stoi(d->bnd));
    2638         735 :   return gsub(gmul(gsqr(gen),ps),gmul2n(gaddgs(gen,d->nbi),2));
    2639             : }
    2640             : 
    2641             : static long
    2642          63 : compute_nbgen(GEN C, long bnd, long nbi)
    2643             : {
    2644             :   struct computeG d;
    2645          63 :   d.C = shifti(C, 1);
    2646          63 :   d.bnd = bnd;
    2647          63 :   d.nbi = nbi;
    2648          63 :   return itos(ground(zbrent((void*)&d, _computeG, gen_2, stoi(bnd), DEFAULTPREC)));
    2649             : }
    2650             : 
    2651             : static GEN
    2652        1714 : _psi(void*E, GEN y)
    2653             : {
    2654        1714 :   GEN lx = (GEN) E;
    2655        1714 :   long prec = realprec(lx);
    2656        1714 :   GEN ly = glog(y, prec);
    2657        1714 :   GEN u = gdiv(lx, ly);
    2658        1714 :   return gsub(gdiv(y ,ly), gpow(u, u, prec));
    2659             : }
    2660             : 
    2661             : static GEN
    2662          63 : opt_param(GEN x, long prec)
    2663             : {
    2664          63 :   return zbrent((void*)glog(x,prec), _psi, gen_2, x, prec);
    2665             : }
    2666             : 
    2667             : static GEN
    2668          63 : check_kernel(long nbg, long N, long prmax, GEN C, GEN M, GEN p, GEN m)
    2669             : {
    2670          63 :   pari_sp av = avma;
    2671          63 :   long lM = lg(M)-1, nbcol = lM;
    2672          63 :   long tbs = maxss(1, expu(nbg/expi(m)));
    2673             :   for (;;)
    2674          14 :   {
    2675          77 :     GEN K = FpMs_leftkernel_elt_col(M, nbcol, N, m);
    2676             :     GEN tab;
    2677          77 :     long i, f=0;
    2678          77 :     long l = lg(K), lm = lgefint(m);
    2679          77 :     GEN idx = diviiexact(subiu(p,1),m), g;
    2680             :     pari_timer ti;
    2681          77 :     if (DEBUGLEVEL) timer_start(&ti);
    2682         154 :     for(i=1; i<l; i++)
    2683         154 :       if (signe(gel(K,i)))
    2684          77 :         break;
    2685          77 :     g = Fp_pow(utoi(i), idx, p);
    2686          77 :     tab = Fp_pow_init(g, p, tbs, p);
    2687          77 :     K = FpC_Fp_mul(K, Fp_inv(gel(K,i), m), m);
    2688      128464 :     for(i=1; i<l; i++)
    2689             :     {
    2690      128387 :       GEN k = gel(K,i);
    2691      128387 :       GEN j = i<=prmax ? utoi(i): addis(C,i-(prmax+1));
    2692      128387 :       if (signe(k)==0 || !equalii(Fp_pow_table(tab, k, p), Fp_pow(j, idx, p)))
    2693       76391 :         gel(K,i) = cgetineg(lm);
    2694             :       else
    2695       51996 :         f++;
    2696             :     }
    2697          77 :     if (DEBUGLEVEL) timer_printf(&ti,"found %ld/%ld logs", f, nbg);
    2698          77 :     if(f > (nbg>>1)) return gerepileupto(av, K);
    2699        4585 :     for(i=1; i<=nbcol; i++)
    2700             :     {
    2701        4571 :       long a = 1+random_Fl(lM);
    2702        4571 :       swap(gel(M,a),gel(M,i));
    2703             :     }
    2704          14 :     if (4*nbcol>5*nbg) nbcol = nbcol*9/10;
    2705             :   }
    2706             : }
    2707             : 
    2708             : static GEN
    2709         126 : Fp_log_find_ind(GEN a, GEN K, long prmax, GEN C, GEN p, GEN m)
    2710             : {
    2711         126 :   pari_sp av=avma;
    2712         126 :   GEN aa = gen_1;
    2713         126 :   long AV = 0;
    2714             :   for(;;)
    2715           0 :   {
    2716         126 :     GEN A = Fp_log_find_rel(a, prmax, C, p, &aa, &AV);
    2717         126 :     GEN F = gel(A,1), E = gel(A,2);
    2718         126 :     GEN Ao = gen_0;
    2719         126 :     long i, l = lg(F);
    2720         962 :     for(i=1; i<l; i++)
    2721             :     {
    2722         836 :       GEN Ki = gel(K,F[i]);
    2723         836 :       if (signe(Ki)<0) break;
    2724         836 :       Ao = addii(Ao, mulis(Ki, E[i]));
    2725             :     }
    2726         126 :     if (i==l) return Fp_divu(Ao, AV, m);
    2727           0 :     aa = gerepileuptoint(av, aa);
    2728             :   }
    2729             : }
    2730             : 
    2731             : static GEN
    2732          63 : Fp_log_index(GEN a, GEN b, GEN m, GEN p)
    2733             : {
    2734          63 :   pari_sp av = avma, av2;
    2735          63 :   long i, j, nbi, nbr = 0, nbrow, nbg;
    2736             :   GEN C, c, Ci, ci, pi, pr, sz, l, Ao, Bo, K, d, p_1;
    2737             :   pari_timer ti;
    2738             :   struct Fp_log_rel r;
    2739          63 :   ulong bnds = itou(roundr_safe(opt_param(sqrti(p),DEFAULTPREC)));
    2740          63 :   ulong bnd = 4*bnds;
    2741          63 :   if (!bnds || cmpii(sqru(bnds),m)>=0) return NULL;
    2742             : 
    2743          63 :   p_1 = subiu(p,1);
    2744          63 :   if (!is_pm1(gcdii(m,diviiexact(p_1,m))))
    2745           0 :     m = diviiexact(p_1, Z_ppo(p_1, m));
    2746          63 :   pr = primes_upto_zv(bnd);
    2747          63 :   nbi = lg(pr)-1;
    2748          63 :   C = sqrtremi(p, &c);
    2749          63 :   av2 = avma;
    2750       12796 :   for (i = 1; i <= nbi; ++i)
    2751             :   {
    2752       12733 :     ulong lp = pr[i];
    2753       26894 :     while (lp <= bnd)
    2754             :     {
    2755       14161 :       nbr++;
    2756       14161 :       lp *= pr[i];
    2757             :     }
    2758             :   }
    2759          63 :   pi = cgetg(nbr+1,t_VECSMALL);
    2760          63 :   Ci = cgetg(nbr+1,t_VECSMALL);
    2761          63 :   ci = cgetg(nbr+1,t_VECSMALL);
    2762          63 :   sz = cgetg(nbr+1,t_VECSMALL);
    2763       12796 :   for (i = 1, j = 1; i <= nbi; ++i)
    2764             :   {
    2765       12733 :     ulong lp = pr[i], sp = expu(2*lp-1);
    2766       26894 :     while (lp <= bnd)
    2767             :     {
    2768       14161 :       pi[j] = lp;
    2769       14161 :       Ci[j] = umodiu(C, lp);
    2770       14161 :       ci[j] = umodiu(c, lp);
    2771       14161 :       sz[j] = sp;
    2772       14161 :       lp *= pr[i];
    2773       14161 :       j++;
    2774             :     }
    2775             :   }
    2776          63 :   r.nbrel = 0;
    2777          63 :   r.nbgen = compute_nbgen(C, bnd, nbi);
    2778          63 :   r.nbmax = 2*(nbi+r.nbgen);
    2779          63 :   r.rel = cgetg(r.nbmax+1,t_VEC);
    2780          63 :   r.prmax = pr[nbi];
    2781          63 :   if (DEBUGLEVEL)
    2782             :   {
    2783           0 :     err_printf("bnd=%lu Size FB=%ld extra gen=%ld \n", bnd, nbi, r.nbgen);
    2784           0 :     timer_start(&ti);
    2785             :   }
    2786          63 :   nbg = Fp_log_sieve(&r, C, c, Ci, ci, pi, sz);
    2787          63 :   nbrow = r.prmax + nbg;
    2788          63 :   if (DEBUGLEVEL)
    2789             :   {
    2790           0 :     err_printf("\n");
    2791           0 :     timer_printf(&ti," %ld relations, %ld generators", r.nbrel, nbi+nbg);
    2792             :   }
    2793          63 :   setlg(r.rel,r.nbrel+1);
    2794          63 :   r.rel = gerepilecopy(av2, r.rel);
    2795          63 :   K = check_kernel(nbi+nbrow-r.prmax, nbrow, r.prmax, C, r.rel, p, m);
    2796          63 :   if (DEBUGLEVEL) timer_start(&ti);
    2797          63 :   Ao = Fp_log_find_ind(a, K, r.prmax, C, p, m);
    2798          63 :   if (DEBUGLEVEL) timer_printf(&ti," log element");
    2799          63 :   Bo = Fp_log_find_ind(b, K, r.prmax, C, p, m);
    2800          63 :   if (DEBUGLEVEL) timer_printf(&ti," log generator");
    2801          63 :   d = gcdii(Ao,Bo);
    2802          63 :   l = Fp_div(diviiexact(Ao, d) ,diviiexact(Bo, d), m);
    2803          63 :   if (!equalii(a,Fp_pow(b,l,p))) pari_err_BUG("Fp_log_index");
    2804          63 :   return gerepileuptoint(av, l);
    2805             : }
    2806             : 
    2807             : static int
    2808     4472902 : Fp_log_use_index(long e, long p)
    2809             : {
    2810     4472902 :   return (e >= 27 && 20*(p+6)<=e*e);
    2811             : }
    2812             : 
    2813             : /* Trivial cases a = 1, -1. Return x s.t. g^x = a or [] if no such x exist */
    2814             : static GEN
    2815     8258051 : Fp_easylog(void *E, GEN a, GEN g, GEN ord)
    2816             : {
    2817     8258051 :   pari_sp av = avma;
    2818     8258051 :   GEN p = (GEN)E;
    2819             :   /* assume a reduced mod p, p not necessarily prime */
    2820     8258051 :   if (equali1(a)) return gen_0;
    2821             :   /* p > 2 */
    2822     5264538 :   if (equalii(subiu(p,1), a))  /* -1 */
    2823             :   {
    2824             :     pari_sp av2;
    2825             :     GEN t;
    2826     1292351 :     ord = get_arith_Z(ord);
    2827     1292351 :     if (mpodd(ord)) { set_avma(av); return cgetg(1, t_VEC); } /* no solution */
    2828     1292337 :     t = shifti(ord,-1); /* only possible solution */
    2829     1292337 :     av2 = avma;
    2830     1292337 :     if (!equalii(Fp_pow(g, t, p), a)) { set_avma(av); return cgetg(1, t_VEC); }
    2831     1292309 :     set_avma(av2); return gerepileuptoint(av, t);
    2832             :   }
    2833     3972187 :   if (typ(ord)==t_INT && BPSW_psp(p) && Fp_log_use_index(expi(ord),expi(p)))
    2834          63 :     return Fp_log_index(a, g, ord, p);
    2835     3972124 :   return gc_NULL(av); /* not easy */
    2836             : }
    2837             : 
    2838             : GEN
    2839     3803343 : Fp_log(GEN a, GEN g, GEN ord, GEN p)
    2840             : {
    2841     3803343 :   GEN v = get_arith_ZZM(ord);
    2842     3803317 :   GEN F = gmael(v,2,1);
    2843     3803317 :   long lF = lg(F)-1, lmax;
    2844     3803317 :   if (lF == 0) return equali1(a)? gen_0: cgetg(1, t_VEC);
    2845     3803289 :   lmax = expi(gel(F,lF));
    2846     3803291 :   if (BPSW_psp(p) && Fp_log_use_index(lmax,expi(p)))
    2847          91 :     v = mkvec2(gel(v,1),ZM_famat_limit(gel(v,2),int2n(27)));
    2848     3803290 :   return gen_PH_log(a,g,v,(void*)p,&Fp_star);
    2849             : }
    2850             : 
    2851             : static ulong
    2852      125870 : Fl_log_naive(ulong a, ulong g, ulong ord, ulong p)
    2853             : {
    2854      125870 :   ulong i, h=1;
    2855      349360 :   for(i=0; i<ord; i++, h = Fl_mul(h, g, p))
    2856      349361 :     if(a==h) return i;
    2857           0 :   return ~0UL;
    2858             : }
    2859             : 
    2860             : static ulong
    2861       20260 : Fl_log_naive_pre(ulong a, ulong g, ulong ord, ulong p, ulong pi)
    2862             : {
    2863       20260 :   ulong i, h=1;
    2864       51007 :   for(i=0; i<ord; i++, h = Fl_mul_pre(h, g, p, pi))
    2865       51007 :     if(a==h) return i;
    2866           0 :   return ~0UL;
    2867             : }
    2868             : 
    2869             : static ulong
    2870           0 : Fl_log_Fp(ulong a, ulong g, ulong ord, ulong p)
    2871             : {
    2872           0 :   pari_sp av = avma;
    2873           0 :   GEN r = Fp_log(utoi(a),utoi(g),utoi(ord),utoi(p));
    2874           0 :   return gc_ulong(av, typ(r)==t_INT ? itou(r): ~0UL);
    2875             : }
    2876             : 
    2877             : ulong
    2878       20260 : Fl_log_pre(ulong a, ulong g, ulong ord, ulong p, ulong pi)
    2879             : {
    2880       20260 :   if (ord <= 200) return Fl_log_naive_pre(a, g, ord, p, pi);
    2881           0 :   return Fl_log_Fp(a, g, ord, p);
    2882             : }
    2883             : 
    2884             : ulong
    2885      125870 : Fl_log(ulong a, ulong g, ulong ord, ulong p)
    2886             : {
    2887      125870 :   if (ord <= 200)
    2888           0 :   return (p&HIGHMASK) ? Fl_log_naive_pre(a, g, ord, p, get_Fl_red(p))
    2889      125870 :                       : Fl_log_naive(a, g, ord, p);
    2890           0 :   return Fl_log_Fp(a, g, ord, p);
    2891             : }
    2892             : 
    2893             : /* find x such that h = g^x mod N > 1, N = prod_{i <= l} P[i]^E[i], P[i] prime.
    2894             :  * PHI[l] = eulerphi(N / P[l]^E[l]).   Destroys P/E */
    2895             : static GEN
    2896         126 : znlog_rec(GEN h, GEN g, GEN N, GEN P, GEN E, GEN PHI)
    2897             : {
    2898         126 :   long l = lg(P) - 1, e = E[l];
    2899         126 :   GEN p = gel(P, l), phi = gel(PHI,l), pe = e == 1? p: powiu(p, e);
    2900             :   GEN a,b, hp,gp, hpe,gpe, ogpe; /* = order(g mod p^e) | p^(e-1)(p-1) */
    2901             : 
    2902         126 :   if (l == 1) {
    2903          98 :     hpe = h;
    2904          98 :     gpe = g;
    2905             :   } else {
    2906          28 :     hpe = modii(h, pe);
    2907          28 :     gpe = modii(g, pe);
    2908             :   }
    2909         126 :   if (e == 1) {
    2910          42 :     hp = hpe;
    2911          42 :     gp = gpe;
    2912             :   } else {
    2913          84 :     hp = remii(hpe, p);
    2914          84 :     gp = remii(gpe, p);
    2915             :   }
    2916         126 :   if (hp == gen_0 || gp == gen_0) return NULL;
    2917         105 :   if (absequaliu(p, 2))
    2918             :   {
    2919          35 :     GEN n = int2n(e);
    2920          35 :     ogpe = Zp_order(gpe, gen_2, e, n);
    2921          35 :     a = Fp_log(hpe, gpe, ogpe, n);
    2922          35 :     if (typ(a) != t_INT) return NULL;
    2923             :   }
    2924             :   else
    2925             :   { /* Avoid black box groups: (Z/p^2)^* / (Z/p)^* ~ (Z/pZ, +), where DL
    2926             :        is trivial */
    2927             :     /* [order(gp), factor(order(gp))] */
    2928          70 :     GEN v = Fp_factored_order(gp, subiu(p,1), p);
    2929          70 :     GEN ogp = gel(v,1);
    2930          70 :     if (!equali1(Fp_pow(hp, ogp, p))) return NULL;
    2931          70 :     a = Fp_log(hp, gp, v, p);
    2932          70 :     if (typ(a) != t_INT) return NULL;
    2933          70 :     if (e == 1) ogpe = ogp;
    2934             :     else
    2935             :     { /* find a s.t. g^a = h (mod p^e), p odd prime, e > 0, (h,p) = 1 */
    2936             :       /* use p-adic log: O(log p + e) mul*/
    2937             :       long vpogpe, vpohpe;
    2938             : 
    2939          28 :       hpe = Fp_mul(hpe, Fp_pow(gpe, negi(a), pe), pe);
    2940          28 :       gpe = Fp_pow(gpe, ogp, pe);
    2941             :       /* g,h = 1 mod p; compute b s.t. h = g^b */
    2942             : 
    2943             :       /* v_p(order g mod pe) */
    2944          28 :       vpogpe = equali1(gpe)? 0: e - Z_pval(subiu(gpe,1), p);
    2945             :       /* v_p(order h mod pe) */
    2946          28 :       vpohpe = equali1(hpe)? 0: e - Z_pval(subiu(hpe,1), p);
    2947          28 :       if (vpohpe > vpogpe) return NULL;
    2948             : 
    2949          28 :       ogpe = mulii(ogp, powiu(p, vpogpe)); /* order g mod p^e */
    2950          28 :       if (is_pm1(gpe)) return is_pm1(hpe)? a: NULL;
    2951          28 :       b = gdiv(Qp_log(cvtop(hpe, p, e)), Qp_log(cvtop(gpe, p, e)));
    2952          28 :       a = addii(a, mulii(ogp, padic_to_Q(b)));
    2953             :     }
    2954             :   }
    2955             :   /* gp^a = hp => x = a mod ogpe => generalized Pohlig-Hellman strategy */
    2956          91 :   if (l == 1) return a;
    2957             : 
    2958          28 :   N = diviiexact(N, pe); /* make N coprime to p */
    2959          28 :   h = Fp_mul(h, Fp_pow(g, modii(negi(a), phi), N), N);
    2960          28 :   g = Fp_pow(g, modii(ogpe, phi), N);
    2961          28 :   setlg(P, l); /* remove last element */
    2962          28 :   setlg(E, l);
    2963          28 :   b = znlog_rec(h, g, N, P, E, PHI);
    2964          28 :   if (!b) return NULL;
    2965          28 :   return addmulii(a, b, ogpe);
    2966             : }
    2967             : 
    2968             : static GEN
    2969          98 : get_PHI(GEN P, GEN E)
    2970             : {
    2971          98 :   long i, l = lg(P);
    2972          98 :   GEN PHI = cgetg(l, t_VEC);
    2973          98 :   gel(PHI,1) = gen_1;
    2974         126 :   for (i=1; i<l-1; i++)
    2975             :   {
    2976          28 :     GEN t, p = gel(P,i);
    2977          28 :     long e = E[i];
    2978          28 :     t = mulii(powiu(p, e-1), subiu(p,1));
    2979          28 :     if (i > 1) t = mulii(t, gel(PHI,i));
    2980          28 :     gel(PHI,i+1) = t;
    2981             :   }
    2982          98 :   return PHI;
    2983             : }
    2984             : 
    2985             : GEN
    2986         238 : znlog(GEN h, GEN g, GEN o)
    2987             : {
    2988         238 :   pari_sp av = avma;
    2989             :   GEN N, fa, P, E, x;
    2990         238 :   switch (typ(g))
    2991             :   {
    2992          28 :     case t_PADIC:
    2993             :     {
    2994          28 :       GEN p = gel(g,2);
    2995          28 :       long v = valp(g);
    2996          28 :       if (v < 0) pari_err_DIM("znlog");
    2997          28 :       if (v > 0) {
    2998           0 :         long k = gvaluation(h, p);
    2999           0 :         if (k % v) return cgetg(1,t_VEC);
    3000           0 :         k /= v;
    3001           0 :         if (!gequal(h, gpowgs(g,k))) { set_avma(av); return cgetg(1,t_VEC); }
    3002           0 :         set_avma(av); return stoi(k);
    3003             :       }
    3004          28 :       N = gel(g,3);
    3005          28 :       g = Rg_to_Fp(g, N);
    3006          28 :       break;
    3007             :     }
    3008         203 :     case t_INTMOD:
    3009         203 :       N = gel(g,1);
    3010         203 :       g = gel(g,2); break;
    3011           7 :     default: pari_err_TYPE("znlog", g);
    3012             :       return NULL; /* LCOV_EXCL_LINE */
    3013             :   }
    3014         231 :   if (equali1(N)) { set_avma(av); return gen_0; }
    3015         231 :   h = Rg_to_Fp(h, N);
    3016         224 :   if (o) return gerepileupto(av, Fp_log(h, g, o, N));
    3017          98 :   fa = Z_factor(N);
    3018          98 :   P = gel(fa,1);
    3019          98 :   E = vec_to_vecsmall(gel(fa,2));
    3020          98 :   x = znlog_rec(h, g, N, P, E, get_PHI(P,E));
    3021          98 :   if (!x) { set_avma(av); return cgetg(1,t_VEC); }
    3022          63 :   return gerepileuptoint(av, x);
    3023             : }
    3024             : 
    3025             : GEN
    3026       62790 : Fp_sqrtn(GEN a, GEN n, GEN p, GEN *zeta)
    3027             : {
    3028       62790 :   if (lgefint(p)==3)
    3029             :   {
    3030       62370 :     long nn = itos_or_0(n);
    3031       62370 :     if (nn)
    3032             :     {
    3033       62370 :       ulong pp = p[2];
    3034             :       ulong uz;
    3035       62370 :       ulong r = Fl_sqrtn(umodiu(a,pp),nn,pp, zeta ? &uz:NULL);
    3036       62349 :       if (r==ULONG_MAX) return NULL;
    3037       62293 :       if (zeta) *zeta = utoi(uz);
    3038       62293 :       return utoi(r);
    3039             :     }
    3040             :   }
    3041         420 :   a = modii(a,p);
    3042         420 :   if (!signe(a))
    3043             :   {
    3044           0 :     if (zeta) *zeta = gen_1;
    3045           0 :     if (signe(n) < 0) pari_err_INV("Fp_sqrtn", mkintmod(gen_0,p));
    3046           0 :     return gen_0;
    3047             :   }
    3048         420 :   if (absequaliu(n,2))
    3049             :   {
    3050         224 :     if (zeta) *zeta = subiu(p,1);
    3051         224 :     return signe(n) > 0 ? Fp_sqrt(a,p): Fp_sqrt(Fp_inv(a, p),p);
    3052             :   }
    3053         196 :   return gen_Shanks_sqrtn(a,n,subiu(p,1),zeta,(void*)p,&Fp_star);
    3054             : }
    3055             : 
    3056             : /*********************************************************************/
    3057             : /**                              FACTORIAL                          **/
    3058             : /*********************************************************************/
    3059             : GEN
    3060       75484 : mulu_interval_step(ulong a, ulong b, ulong step)
    3061             : {
    3062       75484 :   pari_sp av = avma;
    3063             :   ulong k, l, N, n;
    3064             :   long lx;
    3065             :   GEN x;
    3066             : 
    3067       75484 :   if (!a) return gen_0;
    3068       75484 :   if (step == 1) return mulu_interval(a, b);
    3069       75484 :   n = 1 + (b-a) / step;
    3070       75484 :   b -= (b-a) % step;
    3071       75484 :   if (n < 61)
    3072             :   {
    3073       74035 :     if (n == 1) return utoipos(a);
    3074       57645 :     x = muluu(a,a+step); if (n == 2) return x;
    3075      497071 :     for (k=a+2*step; k<=b; k+=step) x = mului(k,x);
    3076       46029 :     return gerepileuptoint(av, x);
    3077             :   }
    3078             :   /* step | b-a */
    3079        1449 :   lx = 1; x = cgetg(2 + n/2, t_VEC);
    3080        1449 :   N = b + a;
    3081        1449 :   for (k = a;; k += step)
    3082             :   {
    3083      230195 :     l = N - k; if (l <= k) break;
    3084      228746 :     gel(x,lx++) = muluu(k,l);
    3085             :   }
    3086        1449 :   if (l == k) gel(x,lx++) = utoipos(k);
    3087        1449 :   setlg(x, lx);
    3088        1449 :   return gerepileuptoint(av, ZV_prod(x));
    3089             : }
    3090             : /* return a * (a+1) * ... * b. Assume a <= b  [ note: factoring out powers of 2
    3091             :  * first is slower ... ] */
    3092             : GEN
    3093      182397 : mulu_interval(ulong a, ulong b)
    3094             : {
    3095      182397 :   pari_sp av = avma;
    3096             :   ulong k, l, N, n;
    3097             :   long lx;
    3098             :   GEN x;
    3099             : 
    3100      182397 :   if (!a) return gen_0;
    3101      182397 :   n = b - a + 1;
    3102      182397 :   if (n < 61)
    3103             :   {
    3104      182383 :     if (n == 1) return utoipos(a);
    3105      126103 :     x = muluu(a,a+1); if (n == 2) return x;
    3106      371697 :     for (k=a+2; k<b; k++) x = mului(k,x);
    3107             :     /* avoid k <= b: broken if b = ULONG_MAX */
    3108       92111 :     return gerepileuptoint(av, mului(b,x));
    3109             :   }
    3110          14 :   lx = 1; x = cgetg(2 + n/2, t_VEC);
    3111          14 :   N = b + a;
    3112          14 :   for (k = a;; k++)
    3113             :   {
    3114        7007 :     l = N - k; if (l <= k) break;
    3115        6993 :     gel(x,lx++) = muluu(k,l);
    3116             :   }
    3117          14 :   if (l == k) gel(x,lx++) = utoipos(k);
    3118          14 :   setlg(x, lx);
    3119          14 :   return gerepileuptoint(av, ZV_prod(x));
    3120             : }
    3121             : GEN
    3122         588 : muls_interval(long a, long b)
    3123             : {
    3124         588 :   pari_sp av = avma;
    3125         588 :   long lx, k, l, N, n = b - a + 1;
    3126             :   GEN x;
    3127             : 
    3128         588 :   if (a <= 0 && b >= 0) return gen_0;
    3129         315 :   if (n < 61)
    3130             :   {
    3131         315 :     x = stoi(a);
    3132         553 :     for (k=a+1; k<=b; k++) x = mulsi(k,x);
    3133         315 :     return gerepileuptoint(av, x);
    3134             :   }
    3135           0 :   lx = 1; x = cgetg(2 + n/2, t_VEC);
    3136           0 :   N = b + a;
    3137           0 :   for (k = a;; k++)
    3138             :   {
    3139           0 :     l = N - k; if (l <= k) break;
    3140           0 :     gel(x,lx++) = mulss(k,l);
    3141             :   }
    3142           0 :   if (l == k) gel(x,lx++) = stoi(k);
    3143           0 :   setlg(x, lx);
    3144           0 :   return gerepileuptoint(av, ZV_prod(x));
    3145             : }
    3146             : 
    3147             : GEN
    3148      495888 : mpfact(long n)
    3149             : {
    3150      495888 :   pari_sp av = avma;
    3151             :   GEN a, v;
    3152             :   long k;
    3153      495888 :   if (n <= 12) switch(n)
    3154             :   {
    3155      440002 :     case 0: case 1: return gen_1;
    3156       34554 :     case 2: return gen_2;
    3157        1304 :     case 3: return utoipos(6);
    3158        2133 :     case 4: return utoipos(24);
    3159         780 :     case 5: return utoipos(120);
    3160         486 :     case 6: return utoipos(720);
    3161         387 :     case 7: return utoipos(5040);
    3162         388 :     case 8: return utoipos(40320);
    3163         416 :     case 9: return utoipos(362880);
    3164         661 :     case 10:return utoipos(3628800);
    3165         249 :     case 11:return utoipos(39916800);
    3166         236 :     case 12:return utoipos(479001600);
    3167           0 :     default: pari_err_DOMAIN("factorial", "argument","<",gen_0,stoi(n));
    3168             :   }
    3169       14292 :   v = cgetg(expu(n) + 2, t_VEC);
    3170       14292 :   for (k = 1;; k++)
    3171       71697 :   {
    3172       85989 :     long m = n >> (k-1), l;
    3173       85989 :     if (m <= 2) break;
    3174       71697 :     l = (1 + (n >> k)) | 1;
    3175             :     /* product of odd numbers in ]n / 2^k, 2 / 2^(k-1)] */
    3176       71697 :     a = mulu_interval_step(l, m, 2);
    3177       71696 :     gel(v,k) = k == 1? a: powiu(a, k);
    3178             :   }
    3179       71697 :   a = gel(v,--k); while (--k) a = mulii(a, gel(v,k));
    3180       14292 :   a = shifti(a, factorial_lval(n, 2));
    3181       14292 :   return gerepileuptoint(av, a);
    3182             : }
    3183             : 
    3184             : ulong
    3185       49897 : factorial_Fl(long n, ulong p)
    3186             : {
    3187             :   long k;
    3188             :   ulong v;
    3189       49897 :   if (p <= (ulong)n) return 0;
    3190       49897 :   v = Fl_powu(2, factorial_lval(n, 2), p);
    3191       49966 :   for (k = 1;; k++)
    3192      128676 :   {
    3193      178642 :     long m = n >> (k-1), l, i;
    3194      178642 :     ulong a = 1;
    3195      178642 :     if (m <= 2) break;
    3196      128712 :     l = (1 + (n >> k)) | 1;
    3197             :     /* product of odd numbers in ]n / 2^k, 2 / 2^(k-1)] */
    3198      748448 :     for (i=l; i<=m; i+=2)
    3199      619761 :       a = Fl_mul(a, i, p);
    3200      128687 :     v = Fl_mul(v, k == 1? a: Fl_powu(a, k, p), p);
    3201             :   }
    3202       49930 :   return v;
    3203             : }
    3204             : 
    3205             : GEN
    3206          60 : factorial_Fp(long n, GEN p)
    3207             : {
    3208          60 :   pari_sp av = avma;
    3209             :   long k;
    3210          60 :   GEN v = Fp_powu(gen_2, factorial_lval(n, 2), p);
    3211          60 :   for (k = 1;; k++)
    3212         134 :   {
    3213         194 :     long m = n >> (k-1), l, i;
    3214         194 :     GEN a = gen_1;
    3215         194 :     if (m <= 2) break;
    3216         134 :     l = (1 + (n >> k)) | 1;
    3217             :     /* product of odd numbers in ]n / 2^k, 2 / 2^(k-1)] */
    3218         402 :     for (i=l; i<=m; i+=2)
    3219         268 :       a = Fp_mulu(a, i, p);
    3220         134 :     v = Fp_mul(v, k == 1? a: Fp_powu(a, k, p), p);
    3221         134 :     v = gerepileuptoint(av, v);
    3222             :   }
    3223          60 :   return v;
    3224             : }
    3225             : 
    3226             : /*******************************************************************/
    3227             : /**                      LUCAS & FIBONACCI                        **/
    3228             : /*******************************************************************/
    3229             : static void
    3230          56 : lucas(ulong n, GEN *a, GEN *b)
    3231             : {
    3232             :   GEN z, t, zt;
    3233          56 :   if (!n) { *a = gen_2; *b = gen_1; return; }
    3234          49 :   lucas(n >> 1, &z, &t); zt = mulii(z, t);
    3235          49 :   switch(n & 3) {
    3236          14 :     case  0: *a = subiu(sqri(z),2); *b = subiu(zt,1); break;
    3237          14 :     case  1: *a = subiu(zt,1);      *b = addiu(sqri(t),2); break;
    3238           7 :     case  2: *a = addiu(sqri(z),2); *b = addiu(zt,1); break;
    3239          14 :     case  3: *a = addiu(zt,1);      *b = subiu(sqri(t),2);
    3240             :   }
    3241          49 : }
    3242             : 
    3243             : GEN
    3244           7 : fibo(long n)
    3245             : {
    3246           7 :   pari_sp av = avma;
    3247             :   GEN a, b;
    3248           7 :   if (!n) return gen_0;
    3249           7 :   lucas((ulong)(labs(n)-1), &a, &b);
    3250           7 :   a = diviuexact(addii(shifti(a,1),b), 5);
    3251           7 :   if (n < 0 && !odd(n)) setsigne(a, -1);
    3252           7 :   return gerepileuptoint(av, a);
    3253             : }
    3254             : 
    3255             : /*******************************************************************/
    3256             : /*                      CONTINUED FRACTIONS                        */
    3257             : /*******************************************************************/
    3258             : static GEN
    3259     2830727 : icopy_lg(GEN x, long l)
    3260             : {
    3261     2830727 :   long lx = lgefint(x);
    3262             :   GEN y;
    3263             : 
    3264     2830727 :   if (lx >= l) return icopy(x);
    3265          35 :   y = cgeti(l); affii(x, y); return y;
    3266             : }
    3267             : 
    3268             : /* continued fraction of a/b. If y != NULL, stop when partial quotients
    3269             :  * differ from y */
    3270             : static GEN
    3271     2831060 : Qsfcont(GEN a, GEN b, GEN y, ulong k)
    3272             : {
    3273             :   GEN  z, c;
    3274     2831060 :   ulong i, l, ly = lgefint(b);
    3275             : 
    3276             :   /* times 1 / log2( (1+sqrt(5)) / 2 )  */
    3277     2831060 :   l = (ulong)(3 + bit_accuracy_mul(ly, 1.44042009041256));
    3278     2831060 :   if (k > 0 && k+1 > 0 && l > k+1) l = k+1; /* beware overflow */
    3279     2831060 :   if (l > LGBITS) l = LGBITS;
    3280             : 
    3281     2831060 :   z = cgetg(l,t_VEC);
    3282     2831060 :   l--;
    3283     2831060 :   if (y) {
    3284         333 :     pari_sp av = avma;
    3285         333 :     if (l >= (ulong)lg(y)) l = lg(y)-1;
    3286       27537 :     for (i = 1; i <= l; i++)
    3287             :     {
    3288       27327 :       GEN q = gel(y,i);
    3289       27327 :       gel(z,i) = q;
    3290       27327 :       c = b; if (!gequal1(q)) c = mulii(q, b);
    3291       27327 :       c = subii(a, c);
    3292       27327 :       if (signe(c) < 0)
    3293             :       { /* partial quotient too large */
    3294         117 :         c = addii(c, b);
    3295         117 :         if (signe(c) >= 0) i++; /* by 1 */
    3296         117 :         break;
    3297             :       }
    3298       27210 :       if (cmpii(c, b) >= 0)
    3299             :       { /* partial quotient too small */
    3300           6 :         c = subii(c, b);
    3301           6 :         if (cmpii(c, b) < 0) {
    3302             :           /* by 1. If next quotient is 1 in y, add 1 */
    3303           6 :           if (i < l && equali1(gel(y,i+1))) gel(z,i) = addiu(q,1);
    3304           6 :           i++;
    3305             :         }
    3306           6 :         break;
    3307             :       }
    3308       27204 :       if ((i & 0xff) == 0) gerepileall(av, 2, &b, &c);
    3309       27204 :       a = b; b = c;
    3310             :     }
    3311             :   } else {
    3312     2830727 :     a = icopy_lg(a, ly);
    3313     2830727 :     b = icopy(b);
    3314    23450113 :     for (i = 1; i <= l; i++)
    3315             :     {
    3316    23449804 :       gel(z,i) = truedvmdii(a,b,&c);
    3317    23449804 :       if (c == gen_0) { i++; break; }
    3318    20619386 :       affii(c, a); cgiv(c); c = a;
    3319    20619386 :       a = b; b = c;
    3320             :     }
    3321             :   }
    3322     2831060 :   i--;
    3323     2831060 :   if (i > 1 && gequal1(gel(z,i)))
    3324             :   {
    3325         100 :     cgiv(gel(z,i)); --i;
    3326         100 :     gel(z,i) = addui(1, gel(z,i)); /* unclean: leave old z[i] on stack */
    3327             :   }
    3328     2831060 :   setlg(z,i+1); return z;
    3329             : }
    3330             : 
    3331             : static GEN
    3332           0 : sersfcont(GEN a, GEN b, long k)
    3333             : {
    3334           0 :   long i, l = typ(a) == t_POL? lg(a): 3;
    3335             :   GEN y, c;
    3336           0 :   if (lg(b) > l) l = lg(b);
    3337           0 :   if (k > 0 && l > k+1) l = k+1;
    3338           0 :   y = cgetg(l,t_VEC);
    3339           0 :   for (i=1; i<l; i++)
    3340             :   {
    3341           0 :     gel(y,i) = poldivrem(a,b,&c);
    3342           0 :     if (gequal0(c)) { i++; break; }
    3343           0 :     a = b; b = c;
    3344             :   }
    3345           0 :   setlg(y, i); return y;
    3346             : }
    3347             : 
    3348             : GEN
    3349     2831742 : gboundcf(GEN x, long k)
    3350             : {
    3351             :   pari_sp av;
    3352     2831742 :   long tx = typ(x), e;
    3353             :   GEN y, a, b, c;
    3354             : 
    3355     2831742 :   if (k < 0) pari_err_DOMAIN("gboundcf","nmax","<",gen_0,stoi(k));
    3356     2831735 :   if (is_scalar_t(tx))
    3357             :   {
    3358     2831735 :     if (gequal0(x)) return mkvec(gen_0);
    3359     2831630 :     switch(tx)
    3360             :     {
    3361         896 :       case t_INT: return mkveccopy(x);
    3362         340 :       case t_REAL:
    3363         340 :         av = avma;
    3364         340 :         c = mantissa_real(x,&e);
    3365         340 :         if (e < 0) pari_err_PREC("gboundcf");
    3366         333 :         y = int2n(e);
    3367         333 :         a = Qsfcont(c,y, NULL, k);
    3368         333 :         b = addsi(signe(x), c);
    3369         333 :         return gerepilecopy(av, Qsfcont(b,y, a, k));
    3370             : 
    3371     2830394 :       case t_FRAC:
    3372     2830394 :         av = avma;
    3373     2830394 :         return gerepileupto(av, Qsfcont(gel(x,1),gel(x,2), NULL, k));
    3374             :     }
    3375           0 :     pari_err_TYPE("gboundcf",x);
    3376             :   }
    3377             : 
    3378           0 :   switch(tx)
    3379             :   {
    3380           0 :     case t_POL: return mkveccopy(x);
    3381           0 :     case t_SER:
    3382           0 :       av = avma;
    3383           0 :       return gerepileupto(av, gboundcf(ser2rfrac_i(x), k));
    3384           0 :     case t_RFRAC:
    3385           0 :       av = avma;
    3386           0 :       return gerepilecopy(av, sersfcont(gel(x,1), gel(x,2), k));
    3387             :   }
    3388           0 :   pari_err_TYPE("gboundcf",x);
    3389             :   return NULL; /* LCOV_EXCL_LINE */
    3390             : }
    3391             : 
    3392             : static GEN
    3393          14 : sfcont2(GEN b, GEN x, long k)
    3394             : {
    3395          14 :   pari_sp av = avma;
    3396          14 :   long lb = lg(b), tx = typ(x), i;
    3397             :   GEN y,p1;
    3398             : 
    3399          14 :   if (k)
    3400             :   {
    3401           7 :     if (k >= lb) pari_err_DIM("contfrac [too few denominators]");
    3402           0 :     lb = k+1;
    3403             :   }
    3404           7 :   y = cgetg(lb,t_VEC);
    3405           7 :   if (lb==1) return y;
    3406           7 :   if (is_scalar_t(tx))
    3407             :   {
    3408           7 :     if (!is_intreal_t(tx) && tx != t_FRAC) pari_err_TYPE("sfcont2",x);
    3409             :   }
    3410           0 :   else if (tx == t_SER) x = ser2rfrac_i(x);
    3411             : 
    3412           7 :   if (!gequal1(gel(b,1))) x = gmul(gel(b,1),x);
    3413           7 :   for (i = 1;;)
    3414             :   {
    3415          35 :     if (tx == t_REAL)
    3416             :     {
    3417          35 :       long e = expo(x);
    3418          35 :       if (e > 0 && nbits2prec(e+1) > realprec(x)) break;
    3419          35 :       gel(y,i) = floorr(x);
    3420          35 :       p1 = subri(x, gel(y,i));
    3421             :     }
    3422             :     else
    3423             :     {
    3424           0 :       gel(y,i) = gfloor(x);
    3425           0 :       p1 = gsub(x, gel(y,i));
    3426             :     }
    3427          35 :     if (++i >= lb) break;
    3428          28 :     if (gequal0(p1)) break;
    3429          28 :     x = gdiv(gel(b,i),p1);
    3430             :   }
    3431           7 :   setlg(y,i);
    3432           7 :   return gerepilecopy(av,y);
    3433             : }
    3434             : 
    3435             : GEN
    3436         112 : gcf(GEN x) { return gboundcf(x,0); }
    3437             : GEN
    3438           0 : gcf2(GEN b, GEN x) { return contfrac0(x,b,0); }
    3439             : GEN
    3440          49 : contfrac0(GEN x, GEN b, long nmax)
    3441             : {
    3442             :   long tb;
    3443             : 
    3444          49 :   if (!b) return gboundcf(x,nmax);
    3445          28 :   tb = typ(b);
    3446          28 :   if (tb == t_INT) return gboundcf(x,itos(b));
    3447          21 :   if (! is_vec_t(tb)) pari_err_TYPE("contfrac0",b);
    3448          21 :   if (nmax < 0) pari_err_DOMAIN("contfrac","nmax","<",gen_0,stoi(nmax));
    3449          14 :   return sfcont2(b,x,nmax);
    3450             : }
    3451             : 
    3452             : GEN
    3453         252 : contfracpnqn(GEN x, long n)
    3454             : {
    3455         252 :   pari_sp av = avma;
    3456         252 :   long i, lx = lg(x);
    3457             :   GEN M,A,B, p0,p1, q0,q1;
    3458             : 
    3459         252 :   if (lx == 1)
    3460             :   {
    3461          28 :     if (! is_matvec_t(typ(x))) pari_err_TYPE("pnqn",x);
    3462          21 :     if (n >= 0) return cgetg(1,t_MAT);
    3463           7 :     return matid(2);
    3464             :   }
    3465         224 :   switch(typ(x))
    3466             :   {
    3467         182 :     case t_VEC: case t_COL: A = x; B = NULL; break;
    3468          42 :     case t_MAT:
    3469          42 :       switch(lgcols(x))
    3470             :       {
    3471           0 :         case 2: A = row(x,1); B = NULL; break;
    3472          35 :         case 3: A = row(x,2); B = row(x,1); break;
    3473           7 :         default: pari_err_DIM("pnqn [ nbrows != 1,2 ]");
    3474             :                  return NULL; /*LCOV_EXCL_LINE*/
    3475             :       }
    3476          35 :       break;
    3477           0 :     default: pari_err_TYPE("pnqn",x);
    3478             :       return NULL; /*LCOV_EXCL_LINE*/
    3479             :   }
    3480         217 :   p1 = gel(A,1);
    3481         217 :   q1 = B? gel(B,1): gen_1; /* p[0], q[0] */
    3482         217 :   if (n >= 0)
    3483             :   {
    3484         182 :     lx = minss(lx, n+2);
    3485         182 :     if (lx == 2) return gerepilecopy(av, mkmat(mkcol2(p1,q1)));
    3486             :   }
    3487          35 :   else if (lx == 2)
    3488           7 :     return gerepilecopy(av, mkmat2(mkcol2(p1,q1), mkcol2(gen_1,gen_0)));
    3489             :   /* lx >= 3 */
    3490         119 :   p0 = gen_1;
    3491         119 :   q0 = gen_0; /* p[-1], q[-1] */
    3492         119 :   M = cgetg(lx, t_MAT);
    3493         119 :   gel(M,1) = mkcol2(p1,q1);
    3494         399 :   for (i=2; i<lx; i++)
    3495             :   {
    3496         280 :     GEN a = gel(A,i), p2,q2;
    3497         280 :     if (B) {
    3498          84 :       GEN b = gel(B,i);
    3499          84 :       p0 = gmul(b,p0);
    3500          84 :       q0 = gmul(b,q0);
    3501             :     }
    3502         280 :     p2 = gadd(gmul(a,p1),p0); p0=p1; p1=p2;
    3503         280 :     q2 = gadd(gmul(a,q1),q0); q0=q1; q1=q2;
    3504         280 :     gel(M,i) = mkcol2(p1,q1);
    3505             :   }
    3506         119 :   if (n < 0) M = mkmat2(gel(M,lx-1), gel(M,lx-2));
    3507         119 :   return gerepilecopy(av, M);
    3508             : }
    3509             : GEN
    3510           0 : pnqn(GEN x) { return contfracpnqn(x,-1); }
    3511             : /* x = [a0, ..., an] from gboundcf, n >= 0;
    3512             :  * return [[p0, ..., pn], [q0,...,qn]] */
    3513             : GEN
    3514      609308 : ZV_allpnqn(GEN x)
    3515             : {
    3516      609308 :   long i, lx = lg(x);
    3517      609308 :   GEN p0, p1, q0, q1, p2, q2, P,Q, v = cgetg(3,t_VEC);
    3518             : 
    3519      609308 :   gel(v,1) = P = cgetg(lx, t_VEC);
    3520      609308 :   gel(v,2) = Q = cgetg(lx, t_VEC);
    3521      609308 :   p0 = gen_1; q0 = gen_0;
    3522      609308 :   gel(P, 1) = p1 = gel(x,1); gel(Q, 1) = q1 = gen_1;
    3523     2092209 :   for (i=2; i<lx; i++)
    3524             :   {
    3525     1482901 :     GEN a = gel(x,i);
    3526     1482901 :     gel(P, i) = p2 = addmulii(p0, a, p1); p0 = p1; p1 = p2;
    3527     1482901 :     gel(Q, i) = q2 = addmulii(q0, a, q1); q0 = q1; q1 = q2;
    3528             :   }
    3529      609308 :   return v;
    3530             : }
    3531             : 
    3532             : /* write Mod(x,N) as a/b, gcd(a,b) = 1, b <= B (no condition if B = NULL) */
    3533             : static GEN
    3534          42 : mod_to_frac(GEN x, GEN N, GEN B)
    3535             : {
    3536             :   GEN a, b, A;
    3537          42 :   if (B) A = divii(shifti(N, -1), B);
    3538             :   else
    3539             :   {
    3540          14 :     A = sqrti(shifti(N, -1));
    3541          14 :     B = A;
    3542             :   }
    3543          42 :   if (!Fp_ratlift(x, N, A,B,&a,&b) || !equali1( gcdii(a,b) )) return NULL;
    3544          28 :   return equali1(b)? a: mkfrac(a,b);
    3545             : }
    3546             : 
    3547             : static GEN
    3548          70 : mod_to_rfrac(GEN x, GEN N, long B)
    3549             : {
    3550             :   GEN a, b;
    3551          70 :   long A, d = degpol(N);
    3552          70 :   if (B >= 0) A = d-1 - B;
    3553             :   else
    3554             :   {
    3555          42 :     B = d >> 1;
    3556          42 :     A = odd(d)? B : B-1;
    3557             :   }
    3558          70 :   if (varn(N) != varn(x)) x = scalarpol(x, varn(N));
    3559          70 :   if (!RgXQ_ratlift(x, N, A, B, &a,&b) || degpol(RgX_gcd(a,b)) > 0) return NULL;
    3560          56 :   return gdiv(a,b);
    3561             : }
    3562             : 
    3563             : /* k > 0 t_INT, x a t_FRAC, returns the convergent a/b
    3564             :  * of the continued fraction of x with b <= k maximal */
    3565             : static GEN
    3566           7 : bestappr_frac(GEN x, GEN k)
    3567             : {
    3568             :   pari_sp av;
    3569             :   GEN p0, p1, p, q0, q1, q, a, y;
    3570             : 
    3571           7 :   if (cmpii(gel(x,2),k) <= 0) return gcopy(x);
    3572           0 :   av = avma; y = x;
    3573           0 :   p1 = gen_1; p0 = truedvmdii(gel(x,1), gel(x,2), &a); /* = floor(x) */
    3574           0 :   q1 = gen_0; q0 = gen_1;
    3575           0 :   x = mkfrac(a, gel(x,2)); /* = frac(x); now 0<= x < 1 */
    3576             :   for(;;)
    3577             :   {
    3578           0 :     x = ginv(x); /* > 1 */
    3579           0 :     a = typ(x)==t_INT? x: divii(gel(x,1), gel(x,2));
    3580           0 :     if (cmpii(a,k) > 0)
    3581             :     { /* next partial quotient will overflow limits */
    3582             :       GEN n, d;
    3583           0 :       a = divii(subii(k, q1), q0);
    3584           0 :       p = addii(mulii(a,p0), p1); p1=p0; p0=p;
    3585           0 :       q = addii(mulii(a,q0), q1); q1=q0; q0=q;
    3586             :       /* compare |y-p0/q0|, |y-p1/q1| */
    3587           0 :       n = gel(y,1);
    3588           0 :       d = gel(y,2);
    3589           0 :       if (abscmpii(mulii(q1, subii(mulii(q0,n), mulii(d,p0))),
    3590             :                    mulii(q0, subii(mulii(q1,n), mulii(d,p1)))) < 0)
    3591           0 :                    { p1 = p0; q1 = q0; }
    3592           0 :       break;
    3593             :     }
    3594           0 :     p = addii(mulii(a,p0), p1); p1=p0; p0=p;
    3595           0 :     q = addii(mulii(a,q0), q1); q1=q0; q0=q;
    3596             : 
    3597           0 :     if (cmpii(q0,k) > 0) break;
    3598           0 :     x = gsub(x,a); /* 0 <= x < 1 */
    3599           0 :     if (typ(x) == t_INT) { p1 = p0; q1 = q0; break; } /* x = 0 */
    3600             : 
    3601             :   }
    3602           0 :   return gerepileupto(av, gdiv(p1,q1));
    3603             : }
    3604             : /* k > 0 t_INT, x != 0 a t_REAL, returns the convergent a/b
    3605             :  * of the continued fraction of x with b <= k maximal */
    3606             : static GEN
    3607     1428430 : bestappr_real(GEN x, GEN k)
    3608             : {
    3609     1428430 :   pari_sp av = avma;
    3610     1428430 :   GEN kr, p0, p1, p, q0, q1, q, a, y = x;
    3611             : 
    3612     1428430 :   p1 = gen_1; a = p0 = floorr(x);
    3613     1428333 :   q1 = gen_0; q0 = gen_1;
    3614     1428333 :   x = subri(x,a); /* 0 <= x < 1 */
    3615     1428325 :   if (!signe(x)) { cgiv(x); return a; }
    3616     1324976 :   kr = itor(k, realprec(x));
    3617             :   for(;;)
    3618     1471032 :   {
    3619             :     long d;
    3620     2796033 :     x = invr(x); /* > 1 */
    3621     2795915 :     if (cmprr(x,kr) > 0)
    3622             :     { /* next partial quotient will overflow limits */
    3623     1312453 :       a = divii(subii(k, q1), q0);
    3624     1312371 :       p = addii(mulii(a,p0), p1); p1=p0; p0=p;
    3625     1312377 :       q = addii(mulii(a,q0), q1); q1=q0; q0=q;
    3626             :       /* compare |y-p0/q0|, |y-p1/q1| */
    3627     1312314 :       if (abscmprr(mulir(q1, subri(mulir(q0,y), p0)),
    3628             :                    mulir(q0, subri(mulir(q1,y), p1))) < 0)
    3629      124576 :                    { p1 = p0; q1 = q0; }
    3630     1312457 :       break;
    3631             :     }
    3632     1483583 :     d = nbits2prec(expo(x) + 1);
    3633     1483583 :     if (d > lg(x)) { p1 = p0; q1 = q0; break; } /* original x was ~ 0 */
    3634             : 
    3635     1483388 :     a = truncr(x); /* truncr(x) will NOT raise e_PREC */
    3636     1483351 :     p = addii(mulii(a,p0), p1); p1=p0; p0=p;
    3637     1483371 :     q = addii(mulii(a,q0), q1); q1=q0; q0=q;
    3638             : 
    3639     1483373 :     if (cmpii(q0,k) > 0) break;
    3640     1477460 :     x = subri(x,a); /* 0 <= x < 1 */
    3641     1477463 :     if (!signe(x)) { p1 = p0; q1 = q0; break; }
    3642             :   }
    3643     1324999 :   if (signe(q1) < 0) { togglesign_safe(&p1); togglesign_safe(&q1); }
    3644     1324999 :   return gerepilecopy(av, equali1(q1)? p1: mkfrac(p1,q1));
    3645             : }
    3646             : 
    3647             : /* k t_INT or NULL */
    3648             : static GEN
    3649     2377894 : bestappr_Q(GEN x, GEN k)
    3650             : {
    3651     2377894 :   long lx, tx = typ(x), i;
    3652             :   GEN a, y;
    3653             : 
    3654     2377894 :   switch(tx)
    3655             :   {
    3656          77 :     case t_INT: return icopy(x);
    3657           7 :     case t_FRAC: return k? bestappr_frac(x, k): gcopy(x);
    3658     1669001 :     case t_REAL:
    3659     1669001 :       if (!signe(x)) return gen_0;
    3660             :       /* i <= e iff nbits2lg(e+1) > lg(x) iff floorr(x) fails */
    3661     1428415 :       i = bit_prec(x); if (i <= expo(x)) return NULL;
    3662     1428430 :       return bestappr_real(x, k? k: int2n(i));
    3663             : 
    3664          28 :     case t_INTMOD: {
    3665          28 :       pari_sp av = avma;
    3666          28 :       a = mod_to_frac(gel(x,2), gel(x,1), k); if (!a) return NULL;
    3667          21 :       return gerepilecopy(av, a);
    3668             :     }
    3669          14 :     case t_PADIC: {
    3670          14 :       pari_sp av = avma;
    3671          14 :       long v = valp(x);
    3672          14 :       a = mod_to_frac(gel(x,4), gel(x,3), k); if (!a) return NULL;
    3673           7 :       if (v) a = gmul(a, powis(gel(x,2), v));
    3674           7 :       return gerepilecopy(av, a);
    3675             :     }
    3676             : 
    3677         322 :     case t_COMPLEX: {
    3678         322 :       pari_sp av = avma;
    3679         322 :       y = cgetg(3, t_COMPLEX);
    3680         322 :       gel(y,2) = bestappr(gel(x,2), k);
    3681         322 :       gel(y,1) = bestappr(gel(x,1), k);
    3682         322 :       if (gequal0(gel(y,2))) return gerepileupto(av, gel(y,1));
    3683          91 :       return y;
    3684             :     }
    3685           0 :     case t_SER:
    3686           0 :       if (ser_isexactzero(x)) return gcopy(x);
    3687             :       /* fall through */
    3688             :     case t_POLMOD: case t_POL: case t_RFRAC:
    3689             :     case t_VEC: case t_COL: case t_MAT:
    3690      708445 :       y = cgetg_copy(x, &lx);
    3691      708480 :       if (lontyp[tx] == 1) i = 1; else { y[1] = x[1]; i = 2; }
    3692     3020145 :       for (; i<lx; i++)
    3693             :       {
    3694     2311678 :         a = bestappr_Q(gel(x,i),k); if (!a) return NULL;
    3695     2311665 :         gel(y,i) = a;
    3696             :       }
    3697      708467 :       if (tx == t_POL) return normalizepol(y);
    3698      708453 :       if (tx == t_SER) return normalizeser(y);
    3699      708453 :       return y;
    3700             :   }
    3701           0 :   pari_err_TYPE("bestappr_Q",x);
    3702             :   return NULL; /* LCOV_EXCL_LINE */
    3703             : }
    3704             : 
    3705             : static GEN
    3706          56 : bestappr_ser(GEN x, long B)
    3707             : {
    3708          56 :   long dN, v = valp(x), lx = lg(x);
    3709             :   GEN t;
    3710          56 :   x = normalizepol(ser2pol_i(x, lx));
    3711          56 :   dN = lx-2;
    3712          56 :   if (v > 0)
    3713             :   {
    3714          14 :     x = RgX_shift_shallow(x, v);
    3715          14 :     dN += v;
    3716             :   }
    3717          42 :   else if (v < 0)
    3718             :   {
    3719           7 :     if (B >= 0) B = maxss(B+v, 0);
    3720             :   }
    3721          56 :   t = mod_to_rfrac(x, pol_xn(dN, varn(x)), B);
    3722          56 :   if (!t) return NULL;
    3723          42 :   if (v < 0)
    3724             :   {
    3725             :     GEN a, b;
    3726             :     long vx;
    3727           7 :     if (typ(t) == t_POL) return RgX_mulXn(t, v);
    3728             :     /* t_RFRAC */
    3729           7 :     vx = varn(x);
    3730           7 :     a = gel(t,1);
    3731           7 :     b = gel(t,2);
    3732           7 :     v -= RgX_valrem(b, &b);
    3733           7 :     if (typ(a) == t_POL && varn(a) == vx) v += RgX_valrem(a, &a);
    3734           7 :     if (v < 0) b = RgX_shift(b, -v);
    3735           0 :     else if (v > 0) {
    3736           0 :       if (typ(a) != t_POL || varn(a) != vx) a = scalarpol_shallow(a, vx);
    3737           0 :       a = RgX_shift(a, v);
    3738             :     }
    3739           7 :     t = mkrfraccopy(a, b);
    3740             :   }
    3741          42 :   return t;
    3742             : }
    3743             : static GEN bestappr_RgX(GEN x, long B);
    3744             : /* x t_POLMOD, B >= 0 or < 0 [omit condition on B].
    3745             :  * Look for coprime t_POL a,b, deg(b)<=B, such that a/b = x */
    3746             : static GEN
    3747          77 : bestappr_RgX(GEN x, long B)
    3748             : {
    3749          77 :   long i, lx, tx = typ(x);
    3750             :   GEN y, t;
    3751          77 :   switch(tx)
    3752             :   {
    3753           0 :     case t_INT: case t_REAL: case t_INTMOD: case t_FRAC:
    3754             :     case t_COMPLEX: case t_PADIC: case t_QUAD: case t_POL:
    3755           0 :       return gcopy(x);
    3756             : 
    3757          14 :     case t_RFRAC: {
    3758          14 :       pari_sp av = avma;
    3759          14 :       if (B < 0 || degpol(gel(x,2)) <= B) return gcopy(x);
    3760           7 :       x = rfrac_to_ser_i(x, 2*B+1);
    3761           7 :       t = bestappr_ser(x, B); if (!t) return NULL;
    3762           0 :       return gerepileupto(av, t);
    3763             :     }
    3764          14 :     case t_POLMOD: {
    3765          14 :       pari_sp av = avma;
    3766          14 :       t = mod_to_rfrac(gel(x,2), gel(x,1), B); if (!t) return NULL;
    3767          14 :       return gerepileupto(av, t);
    3768             :     }
    3769          49 :     case t_SER: {
    3770          49 :       pari_sp av = avma;
    3771          49 :       t = bestappr_ser(x, B); if (!t) return NULL;
    3772          42 :       return gerepileupto(av, t);
    3773             :     }
    3774             : 
    3775           0 :     case t_VEC: case t_COL: case t_MAT:
    3776           0 :       y = cgetg_copy(x, &lx);
    3777           0 :       if (lontyp[tx] == 1) i = 1; else { y[1] = x[1]; i = 2; }
    3778           0 :       for (; i<lx; i++)
    3779             :       {
    3780           0 :         t = bestappr_RgX(gel(x,i),B); if (!t) return NULL;
    3781           0 :         gel(y,i) = t;
    3782             :       }
    3783           0 :       return y;
    3784             :   }
    3785           0 :   pari_err_TYPE("bestappr_RgX",x);
    3786             :   return NULL; /* LCOV_EXCL_LINE */
    3787             : }
    3788             : 
    3789             : /* allow k = NULL: maximal accuracy */
    3790             : GEN
    3791       66215 : bestappr(GEN x, GEN k)
    3792             : {
    3793       66215 :   pari_sp av = avma;
    3794       66215 :   if (k) { /* replace by floor(k) */
    3795       65928 :     switch(typ(k))
    3796             :     {
    3797        4424 :       case t_INT:
    3798        4424 :         break;
    3799       61504 :       case t_REAL: case t_FRAC:
    3800       61504 :         k = floor_safe(k); /* left on stack for efficiency */
    3801       61507 :         if (!signe(k)) k = gen_1;
    3802       61507 :         break;
    3803           0 :       default:
    3804           0 :         pari_err_TYPE("bestappr [bound type]", k);
    3805           0 :         break;
    3806             :     }
    3807         287 :   }
    3808       66218 :   x = bestappr_Q(x, k);
    3809       66218 :   if (!x) { set_avma(av); return cgetg(1,t_VEC); }
    3810       66203 :   return x;
    3811             : }
    3812             : GEN
    3813          77 : bestapprPade(GEN x, long B)
    3814             : {
    3815          77 :   pari_sp av = avma;
    3816          77 :   GEN t = bestappr_RgX(x, B);
    3817          77 :   if (!t) { set_avma(av); return cgetg(1,t_VEC); }
    3818          63 :   return t;
    3819             : }

Generated by: LCOV version 1.13