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 - trans1.c (source / functions) Hit Total Coverage
Test: PARI/GP v2.12.0 lcov report (development 23694-b3ccec097) Lines: 2038 2098 97.1 %
Date: 2019-03-20 05:44:21 Functions: 151 153 98.7 %
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. It is distributed in the hope that it will be useful, but WITHOUT
       8             : ANY WARRANTY WHATSOEVER.
       9             : 
      10             : Check the License for details. You should have received a copy of it, along
      11             : with the package; see the file 'COPYING'. If not, write to the Free Software
      12             : Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. */
      13             : 
      14             : /********************************************************************/
      15             : /**                                                                **/
      16             : /**                   TRANSCENDENTAL FUNCTIONS                     **/
      17             : /**                                                                **/
      18             : /********************************************************************/
      19             : #include "pari.h"
      20             : #include "paripriv.h"
      21             : 
      22             : #ifdef LONG_IS_64BIT
      23             : static const long SQRTVERYBIGINT = 3037000500L; /* ceil(sqrt(LONG_MAX)) */
      24             : #else
      25             : static const long SQRTVERYBIGINT = 46341L;
      26             : #endif
      27             : 
      28             : static THREAD GEN gcatalan, geuler, glog2, gpi;
      29             : void
      30      123066 : pari_init_floats(void)
      31             : {
      32      123066 :   gcatalan = geuler = gpi = bernzone = glog2 = NULL;
      33      123066 : }
      34             : 
      35             : void
      36      121797 : pari_close_floats(void)
      37             : {
      38      121797 :   guncloneNULL(gcatalan);
      39      121475 :   guncloneNULL(geuler);
      40      121116 :   guncloneNULL(gpi);
      41      121022 :   guncloneNULL(glog2);
      42      120988 :   guncloneNULL_deep(bernzone);
      43      120739 : }
      44             : 
      45             : /********************************************************************/
      46             : /**                   GENERIC BINARY SPLITTING                     **/
      47             : /**                    (Haible, Papanikolaou)                      **/
      48             : /********************************************************************/
      49             : void
      50        7592 : abpq_init(struct abpq *A, long n)
      51             : {
      52        7592 :   A->a = (GEN*)new_chunk(n+1);
      53        7592 :   A->b = (GEN*)new_chunk(n+1);
      54        7592 :   A->p = (GEN*)new_chunk(n+1);
      55        7592 :   A->q = (GEN*)new_chunk(n+1);
      56        7592 : }
      57             : static GEN
      58      707548 : mulii3(GEN a, GEN b, GEN c) { return mulii(mulii(a,b),c); }
      59             : static GEN
      60      209575 : mulii4(GEN a, GEN b, GEN c, GEN d) { return mulii(mulii(a,b),mulii(c,d)); }
      61             : 
      62             : /* T_{n1,n1+1}, given P = p[n1]p[n1+1] */
      63             : static GEN
      64      209575 : T2(struct abpq *A, long n1, GEN P)
      65             : {
      66      209575 :   GEN u1 = mulii4(A->a[n1], A->p[n1], A->b[n1+1], A->q[n1+1]);
      67      209575 :   GEN u2 = mulii3(A->b[n1],A->a[n1+1], P);
      68      209575 :   return addii(u1, u2);
      69             : }
      70             : 
      71             : /* assume n2 > n1. Compute sum_{n1 <= n < n2} a/b(n) p/q(n1)... p/q(n) */
      72             : void
      73      411648 : abpq_sum(struct abpq_res *r, long n1, long n2, struct abpq *A)
      74             : {
      75             :   struct abpq_res L, R;
      76             :   GEN u1, u2;
      77             :   pari_sp av;
      78             :   long n;
      79      411648 :   switch(n2 - n1)
      80             :   {
      81             :     GEN b, p, q;
      82             :     case 1:
      83          45 :       r->P = A->p[n1];
      84          45 :       r->Q = A->q[n1];
      85          45 :       r->B = A->b[n1];
      86          45 :       r->T = mulii(A->a[n1], A->p[n1]);
      87      209665 :       return;
      88             :     case 2:
      89      115658 :       r->P = mulii(A->p[n1], A->p[n1+1]);
      90      115658 :       r->Q = mulii(A->q[n1], A->q[n1+1]);
      91      115658 :       r->B = mulii(A->b[n1], A->b[n1+1]);
      92      115658 :       av = avma;
      93      115658 :       r->T = gerepileuptoint(av, T2(A, n1, r->P));
      94      115658 :       return;
      95             : 
      96             :     case 3:
      97       93917 :       p = mulii(A->p[n1+1], A->p[n1+2]);
      98       93917 :       q = mulii(A->q[n1+1], A->q[n1+2]);
      99       93917 :       b = mulii(A->b[n1+1], A->b[n1+2]);
     100       93917 :       r->P = mulii(A->p[n1], p);
     101       93917 :       r->Q = mulii(A->q[n1], q);
     102       93917 :       r->B = mulii(A->b[n1], b);
     103       93917 :       av = avma;
     104       93917 :       u1 = mulii3(b, q, A->a[n1]);
     105       93917 :       u2 = mulii(A->b[n1], T2(A, n1+1, p));
     106       93917 :       r->T = gerepileuptoint(av, mulii(A->p[n1], addii(u1, u2)));
     107       93917 :       return;
     108             :   }
     109             : 
     110      202028 :   av = avma;
     111      202028 :   n = (n1 + n2) >> 1;
     112      202028 :   abpq_sum(&L, n1, n, A);
     113      202028 :   abpq_sum(&R, n, n2, A);
     114             : 
     115      202028 :   r->P = mulii(L.P, R.P);
     116      202028 :   r->Q = mulii(L.Q, R.Q);
     117      202028 :   r->B = mulii(L.B, R.B);
     118      202028 :   u1 = mulii3(R.B,R.Q,L.T);
     119      202028 :   u2 = mulii3(L.B,L.P,R.T);
     120      202028 :   r->T = addii(u1,u2);
     121      202028 :   set_avma(av);
     122      202028 :   r->P = icopy(r->P);
     123      202028 :   r->Q = icopy(r->Q);
     124      202028 :   r->B = icopy(r->B);
     125      202028 :   r->T = icopy(r->T);
     126             : }
     127             : 
     128             : /********************************************************************/
     129             : /**                                                                **/
     130             : /**                               PI                               **/
     131             : /**                                                                **/
     132             : /********************************************************************/
     133             : /* replace *old clone by c. Protect against SIGINT */
     134             : static void
     135        5321 : swap_clone(GEN *old, GEN c)
     136        5321 : { GEN tmp = *old; *old = c; guncloneNULL(tmp); }
     137             : 
     138             : /*                         ----
     139             :  *  53360 (640320)^(1/2)   \    (6n)! (545140134 n + 13591409)
     140             :  *  -------------------- = /    ------------------------------
     141             :  *        Pi               ----   (n!)^3 (3n)! (-640320)^(3n)
     142             :  *                         n>=0
     143             :  *
     144             :  * Ramanujan's formula + binary splitting */
     145             : static GEN
     146        2212 : pi_ramanujan(long prec)
     147             : {
     148        2212 :   const ulong B = 545140134, A = 13591409, C = 640320;
     149        2212 :   const double alpha2 = 47.11041314; /* 3log(C/12) / log(2) */
     150             :   long n, nmax, prec2;
     151             :   struct abpq_res R;
     152             :   struct abpq S;
     153             :   GEN D, u;
     154             : 
     155        2212 :   nmax = (long)(1 + prec2nbits(prec)/alpha2);
     156             : #ifdef LONG_IS_64BIT
     157        1892 :   D = utoipos(10939058860032000UL); /* C^3/24 */
     158             : #else
     159         320 :   D = uutoi(2546948UL,495419392UL);
     160             : #endif
     161        2212 :   abpq_init(&S, nmax);
     162        2212 :   S.a[0] = utoipos(A);
     163        2212 :   S.b[0] = S.p[0] = S.q[0] = gen_1;
     164       42647 :   for (n = 1; n <= nmax; n++)
     165             :   {
     166       40435 :     S.a[n] = addiu(muluu(B, n), A);
     167       40435 :     S.b[n] = gen_1;
     168       40435 :     S.p[n] = mulis(muluu(6*n-5, 2*n-1), 1-6*n);
     169       40435 :     S.q[n] = mulii(sqru(n), muliu(D,n));
     170             :   }
     171        2212 :   abpq_sum(&R, 0, nmax, &S); prec2 = prec+EXTRAPRECWORD;
     172        2212 :   u = itor(muliu(R.Q,C/12), prec2);
     173        2212 :   return rtor(mulrr(divri(u, R.T), sqrtr_abs(utor(C,prec2))), prec);
     174             : }
     175             : 
     176             : #if 0 /* Much slower than binary splitting at least up to prec = 10^8 */
     177             : /* Gauss - Brent-Salamin AGM iteration */
     178             : static GEN
     179             : pi_brent_salamin(long prec)
     180             : {
     181             :   GEN A, B, C;
     182             :   pari_sp av2;
     183             :   long i, G;
     184             : 
     185             :   G = - prec2nbits(prec);
     186             :   incrprec(prec);
     187             : 
     188             :   A = real2n(-1, prec);
     189             :   B = sqrtr_abs(A); /* = 1/sqrt(2) */
     190             :   setexpo(A, 0);
     191             :   C = real2n(-2, prec); av2 = avma;
     192             :   for (i = 0;; i++)
     193             :   {
     194             :     GEN y, a, b, B_A = subrr(B, A);
     195             :     pari_sp av3 = avma;
     196             :     if (expo(B_A) < G) break;
     197             :     a = addrr(A,B); shiftr_inplace(a, -1);
     198             :     b = mulrr(A,B);
     199             :     affrr(a, A);
     200             :     affrr(sqrtr_abs(b), B); set_avma(av3);
     201             :     y = sqrr(B_A); shiftr_inplace(y, i - 2);
     202             :     affrr(subrr(C, y), C); set_avma(av2);
     203             :   }
     204             :   shiftr_inplace(C, 2);
     205             :   return divrr(sqrr(addrr(A,B)), C);
     206             : }
     207             : #endif
     208             : 
     209             : GEN
     210    14575849 : constpi(long prec)
     211             : {
     212             :   pari_sp av;
     213             :   GEN tmp;
     214    14575849 :   if (gpi && realprec(gpi) >= prec) return gpi;
     215             : 
     216        2212 :   av = avma;
     217        2212 :   tmp = gclone(pi_ramanujan(prec));
     218        2212 :   swap_clone(&gpi,tmp);
     219        2212 :   set_avma(av); return gpi;
     220             : }
     221             : 
     222             : GEN
     223    14575152 : mppi(long prec) { return rtor(constpi(prec), prec); }
     224             : 
     225             : /* Pi * 2^n */
     226             : GEN
     227     5161938 : Pi2n(long n, long prec)
     228             : {
     229     5161938 :   GEN x = mppi(prec); shiftr_inplace(x, n);
     230     5161938 :   return x;
     231             : }
     232             : 
     233             : /* I * Pi * 2^n */
     234             : GEN
     235        9555 : PiI2n(long n, long prec) { retmkcomplex(gen_0, Pi2n(n, prec)); }
     236             : 
     237             : /* 2I * Pi */
     238             : GEN
     239        4361 : PiI2(long prec) { return PiI2n(1, prec); }
     240             : 
     241             : /********************************************************************/
     242             : /**                                                                **/
     243             : /**                       EULER CONSTANT                           **/
     244             : /**                                                                **/
     245             : /********************************************************************/
     246             : 
     247             : GEN
     248       31420 : consteuler(long prec)
     249             : {
     250             :   GEN u,v,a,b,tmpeuler;
     251             :   long l, n1, n, k, x;
     252             :   pari_sp av1, av2;
     253             : 
     254       31420 :   if (geuler && realprec(geuler) >= prec) return geuler;
     255             : 
     256         412 :   av1 = avma; tmpeuler = cgetr_block(prec);
     257             : 
     258         412 :   incrprec(prec);
     259             : 
     260         412 :   l = prec+EXTRAPRECWORD; x = (long) (1 + prec2nbits_mul(l, M_LN2/4));
     261         412 :   a = utor(x,l); u=logr_abs(a); setsigne(u,-1); affrr(u,a);
     262         412 :   b = real_1(l);
     263         412 :   v = real_1(l);
     264         412 :   n = (long)(1+3.591*x); /* z=3.591: z*[ ln(z)-1 ]=1 */
     265         412 :   n1 = minss(n, SQRTVERYBIGINT);
     266         412 :   if (x < SQRTVERYBIGINT)
     267             :   {
     268         412 :     ulong xx = x*x;
     269         412 :     av2 = avma;
     270      142742 :     for (k=1; k<n1; k++)
     271             :     {
     272      142330 :       affrr(divru(mulur(xx,b),k*k), b);
     273      142330 :       affrr(divru(addrr(divru(mulur(xx,a),k),b),k), a);
     274      142330 :       affrr(addrr(u,a), u);
     275      142330 :       affrr(addrr(v,b), v); set_avma(av2);
     276             :     }
     277         824 :     for (   ; k<=n; k++)
     278             :     {
     279         412 :       affrr(divru(divru(mulur(xx,b),k),k), b);
     280         412 :       affrr(divru(addrr(divru(mulur(xx,a),k),b),k), a);
     281         412 :       affrr(addrr(u,a), u);
     282         412 :       affrr(addrr(v,b), v); set_avma(av2);
     283             :     }
     284             :   }
     285             :   else
     286             :   {
     287           0 :     GEN xx = sqru(x);
     288           0 :     av2 = avma;
     289           0 :     for (k=1; k<n1; k++)
     290             :     {
     291           0 :       affrr(divru(mulir(xx,b),k*k), b);
     292           0 :       affrr(divru(addrr(divru(mulir(xx,a),k),b),k), a);
     293           0 :       affrr(addrr(u,a), u);
     294           0 :       affrr(addrr(v,b), v); set_avma(av2);
     295             :     }
     296           0 :     for (   ; k<=n; k++)
     297             :     {
     298           0 :       affrr(divru(divru(mulir(xx,b),k),k), b);
     299           0 :       affrr(divru(addrr(divru(mulir(xx,a),k),b),k), a);
     300           0 :       affrr(addrr(u,a), u);
     301           0 :       affrr(addrr(v,b), v); set_avma(av2);
     302             :     }
     303             :   }
     304         412 :   divrrz(u,v,tmpeuler);
     305         412 :   swap_clone(&geuler,tmpeuler);
     306         412 :   set_avma(av1); return geuler;
     307             : }
     308             : 
     309             : GEN
     310       31420 : mpeuler(long prec) { return rtor(consteuler(prec), prec); }
     311             : 
     312             : /********************************************************************/
     313             : /**                                                                **/
     314             : /**                       CATALAN CONSTANT                         **/
     315             : /**                                                                **/
     316             : /********************************************************************/
     317             : /* 8G = 3\sum_{n>=0} 1/(binomial(2n,n)(2n+1)^2) + Pi log(2+sqrt(3)) */
     318             : static GEN
     319          14 : catalan(long prec)
     320             : {
     321          14 :   long i, nmax = prec2nbits(prec) >> 1;
     322             :   struct abpq_res R;
     323             :   struct abpq A;
     324             :   GEN u, v;
     325          14 :   abpq_init(&A, nmax);
     326          14 :   A.a[0] = A.b[0] = A.p[0] = A.q[0] = gen_1;
     327        6510 :   for (i = 1; i <= nmax; i++)
     328             :   {
     329        6496 :     A.a[i] = gen_1;
     330        6496 :     A.b[i] = utoipos((i<<1)+1);
     331        6496 :     A.p[i] = utoipos(i);
     332        6496 :     A.q[i] = utoipos((i<<2)+2);
     333             :   }
     334          14 :   abpq_sum(&R, 0, nmax, &A);
     335          14 :   u = mulur(3, rdivii(R.T, mulii(R.B,R.Q),prec));
     336          14 :   v = mulrr(mppi(prec), logr_abs(addrs(sqrtr_abs(utor(3,prec)), 2)));
     337          14 :   u = addrr(u, v); shiftr_inplace(u, -3);
     338          14 :   return u;
     339             : }
     340             : 
     341             : GEN
     342          14 : constcatalan(long prec)
     343             : {
     344          14 :   pari_sp av = avma;
     345             :   GEN tmp;
     346          14 :   if (gcatalan && realprec(gcatalan) >= prec) return gcatalan;
     347          14 :   tmp = gclone(catalan(prec));
     348          14 :   swap_clone(&gcatalan,tmp);
     349          14 :   set_avma(av); return gcatalan;
     350             : }
     351             : 
     352             : GEN
     353          14 : mpcatalan(long prec) { return rtor(constcatalan(prec), prec); }
     354             : 
     355             : /********************************************************************/
     356             : /**                                                                **/
     357             : /**          TYPE CONVERSION FOR TRANSCENDENTAL FUNCTIONS          **/
     358             : /**                                                                **/
     359             : /********************************************************************/
     360             : static GEN
     361       90021 : transvec(GEN (*f)(GEN,long), GEN x, long prec)
     362             : {
     363             :   long i, l;
     364       90021 :   GEN y = cgetg_copy(x, &l);
     365       90021 :   for (i=1; i<l; i++) gel(y,i) = f(gel(x,i),prec);
     366       90014 :   return y;
     367             : }
     368             : 
     369             : GEN
     370      656063 : trans_eval(const char *fun, GEN (*f)(GEN,long), GEN x, long prec)
     371             : {
     372      656063 :   pari_sp av = avma;
     373      656063 :   if (prec < 3) pari_err_BUG("trans_eval [prec < 3]");
     374      656063 :   switch(typ(x))
     375             :   {
     376      415381 :     case t_INT:    x = f(itor(x,prec),prec); break;
     377      150605 :     case t_FRAC:   x = f(fractor(x, prec),prec); break;
     378           7 :     case t_QUAD:   x = f(quadtofp(x,prec),prec); break;
     379          14 :     case t_POLMOD: x = transvec(f, polmod_to_embed(x,prec), prec); break;
     380             :     case t_VEC:
     381             :     case t_COL:
     382       90007 :     case t_MAT: return transvec(f, x, prec);
     383          49 :     default: pari_err_TYPE(fun,x); return NULL;
     384             :   }
     385      565986 :   return gerepileupto(av, x);
     386             : }
     387             : 
     388             : /*******************************************************************/
     389             : /*                                                                 */
     390             : /*                            POWERING                             */
     391             : /*                                                                 */
     392             : /*******************************************************************/
     393             : /* x a t_REAL 0, return exp(x) */
     394             : static GEN
     395      212399 : mpexp0(GEN x)
     396             : {
     397      212399 :   long e = expo(x);
     398      212399 :   return e >= 0? real_0_bit(e): real_1_bit(-e);
     399             : }
     400             : static GEN
     401        3815 : powr0(GEN x)
     402        3815 : { return signe(x)? real_1(realprec(x)): mpexp0(x); }
     403             : 
     404             : /* x t_POL or t_SER, return scalarpol(Rg_get_1(x)) */
     405             : static GEN
     406      177165 : scalarpol_get_1(GEN x)
     407             : {
     408      177165 :   GEN y = cgetg(3,t_POL);
     409      177165 :   y[1] = evalvarn(varn(x)) | evalsigne(1);
     410      177165 :   gel(y,2) = Rg_get_1(x); return y;
     411             : }
     412             : /* to be called by the generic function gpowgs(x,s) when s = 0 */
     413             : static GEN
     414     3698058 : gpowg0(GEN x)
     415             : {
     416             :   long lx, i;
     417             :   GEN y;
     418             : 
     419     3698058 :   switch(typ(x))
     420             :   {
     421             :     case t_INT: case t_REAL: case t_FRAC: case t_PADIC:
     422      809936 :       return gen_1;
     423             : 
     424           7 :     case t_QUAD: x++; /*fall through*/
     425             :     case t_COMPLEX: {
     426       21284 :       pari_sp av = avma;
     427       21284 :       GEN a = gpowg0(gel(x,1));
     428       21284 :       GEN b = gpowg0(gel(x,2));
     429       21284 :       if (a == gen_1) return b;
     430          14 :       if (b == gen_1) return a;
     431           7 :       return gerepileupto(av, gmul(a,b));
     432             :     }
     433             :     case t_INTMOD:
     434         119 :       y = cgetg(3,t_INTMOD);
     435         119 :       gel(y,1) = icopy(gel(x,1));
     436         119 :       gel(y,2) = is_pm1(gel(x,1))? gen_0: gen_1;
     437         119 :       return y;
     438             : 
     439        5761 :     case t_FFELT: return FF_1(x);
     440             : 
     441             :     case t_POLMOD:
     442        1071 :       retmkpolmod(scalarpol_get_1(gel(x,1)), gcopy(gel(x,1)));
     443             : 
     444             :     case t_RFRAC:
     445           7 :       return scalarpol_get_1(gel(x,2));
     446             :     case t_POL: case t_SER:
     447      176087 :       return scalarpol_get_1(x);
     448             : 
     449             :     case t_MAT:
     450          84 :       lx=lg(x); if (lx==1) return cgetg(1,t_MAT);
     451          77 :       if (lx != lgcols(x)) pari_err_DIM("gpow");
     452          77 :       y = matid(lx-1);
     453          77 :       for (i=1; i<lx; i++) gcoeff(y,i,i) = gpowg0(gcoeff(x,i,i));
     454          77 :       return y;
     455          14 :     case t_QFR: return qfr_1(x);
     456     2683681 :     case t_QFI: return qfi_1(x);
     457           7 :     case t_VECSMALL: return identity_perm(lg(x) - 1);
     458             :   }
     459           7 :   pari_err_TYPE("gpow",x);
     460             :   return NULL; /* LCOV_EXCL_LINE */
     461             : }
     462             : 
     463             : static GEN
     464    42626143 : _sqr(void *data /* ignored */, GEN x) { (void)data; return gsqr(x); }
     465             : static GEN
     466    17368688 : _mul(void *data /* ignored */, GEN x, GEN y) { (void)data; return gmul(x,y); }
     467             : static GEN
     468      110700 : _one(void *x) { return gpowg0((GEN) x); }
     469             : static GEN
     470     5352845 : _sqri(void *data /* ignored */, GEN x) { (void)data; return sqri(x); }
     471             : static GEN
     472     3661591 : _muli(void *data /* ignored */, GEN x, GEN y) { (void)data; return mulii(x,y); }
     473             : static GEN
     474     9839878 : _sqrr(void *data /* ignored */, GEN x) { (void)data; return sqrr(x); }
     475             : static GEN
     476     1714217 : _mulr(void *data /* ignored */, GEN x, GEN y) { (void)data; return mulrr(x,y); }
     477             : static GEN
     478       14505 : _oner(void *data /* prec */) { return real_1( *(long*) data); }
     479             : 
     480             : /* INTEGER POWERING (a^n for integer a != 0 and integer n > 0)
     481             :  *
     482             :  * Use left shift binary algorithm (RS is wasteful: multiplies big numbers,
     483             :  * with LS one of them is the base, hence small). Sign of result is set
     484             :  * to s (= 1,-1). Makes life easier for caller, which otherwise might do a
     485             :  * setsigne(gen_1 / gen_m1) */
     486             : static GEN
     487    39918155 : powiu_sign(GEN a, ulong N, long s)
     488             : {
     489             :   pari_sp av;
     490             :   GEN y;
     491             : 
     492    39918155 :   if (lgefint(a) == 3)
     493             :   { /* easy if |a| < 3 */
     494    36082128 :     ulong q = a[2];
     495    36082128 :     if (q == 1) return (s>0)? gen_1: gen_m1;
     496    27570229 :     if (q == 2) { a = int2u(N); setsigne(a,s); return a; }
     497    16152191 :     q = upowuu(q, N);
     498    16152290 :     if (q) return s>0? utoipos(q): utoineg(q);
     499             :   }
     500     4512780 :   if (N <= 2) {
     501     2888622 :     if (N == 2) return sqri(a);
     502       23936 :     a = icopy(a); setsigne(a,s); return a;
     503             :   }
     504     1624158 :   av = avma;
     505     1624158 :   y = gen_powu_i(a, N, NULL, &_sqri, &_muli);
     506     1624158 :   setsigne(y,s); return gerepileuptoint(av, y);
     507             : }
     508             : /* a^n */
     509             : GEN
     510    39828086 : powiu(GEN a, ulong n)
     511             : {
     512             :   long s;
     513    39828086 :   if (!n) return gen_1;
     514    39365658 :   s = signe(a);
     515    39365658 :   if (!s) return gen_0;
     516    39352518 :   return powiu_sign(a, n, (s < 0 && odd(n))? -1: 1);
     517             : }
     518             : GEN
     519    22334699 : powis(GEN a, long n)
     520             : {
     521             :   long s;
     522             :   GEN t, y;
     523    22334699 :   if (n >= 0) return powiu(a, n);
     524      109623 :   s = signe(a);
     525      109623 :   if (!s) pari_err_INV("powis",gen_0);
     526      109623 :   t = (s < 0 && odd(n))? gen_m1: gen_1;
     527      109623 :   if (is_pm1(a)) return t;
     528             :   /* n < 0, |a| > 1 */
     529      109329 :   y = cgetg(3,t_FRAC);
     530      109329 :   gel(y,1) = t;
     531      109329 :   gel(y,2) = powiu_sign(a, -n, 1); /* force denominator > 0 */
     532      109329 :   return y;
     533             : }
     534             : GEN
     535    33918085 : powuu(ulong p, ulong N)
     536             : {
     537    33918085 :   pari_sp av = avma;
     538    33918085 :   long P[] = {evaltyp(t_INT)|_evallg(3), evalsigne(1)|evallgefint(3),0};
     539             :   ulong pN;
     540             :   GEN y;
     541    33918085 :   if (N <= 2)
     542             :   {
     543    28649177 :     if (N == 2) return sqru(p);
     544    24842307 :     if (N == 1) return utoi(p);
     545     4861622 :     return gen_1;
     546             :   }
     547     5268908 :   if (!p) return gen_0;
     548     5268908 :   pN = upowuu(p, N);
     549     5268908 :   if (pN) return utoipos(pN);
     550      814922 :   if (p == 2) return int2u(N);
     551      810324 :   P[2] = p; av = avma;
     552      810324 :   y = gen_powu_i(P, N, NULL, &_sqri, &_muli);
     553      810324 :   return gerepileuptoint(av, y);
     554             : }
     555             : 
     556             : /* return 0 if overflow */
     557             : static ulong
     558     5341658 : usqru(ulong p) { return p & HIGHMASK? 0: p*p; }
     559             : ulong
     560    28222811 : upowuu(ulong p, ulong k)
     561             : {
     562             : #ifdef LONG_IS_64BIT
     563    24187341 :   const ulong CUTOFF3 = 2642245;
     564    24187341 :   const ulong CUTOFF4 = 65535;
     565    24187341 :   const ulong CUTOFF5 = 7131;
     566    24187341 :   const ulong CUTOFF6 = 1625;
     567    24187341 :   const ulong CUTOFF7 = 565;
     568    24187341 :   const ulong CUTOFF8 = 255;
     569    24187341 :   const ulong CUTOFF9 = 138;
     570    24187341 :   const ulong CUTOFF10 = 84;
     571    24187341 :   const ulong CUTOFF11 = 56;
     572    24187341 :   const ulong CUTOFF12 = 40;
     573    24187341 :   const ulong CUTOFF13 = 30;
     574    24187341 :   const ulong CUTOFF14 = 23;
     575    24187341 :   const ulong CUTOFF15 = 19;
     576    24187341 :   const ulong CUTOFF16 = 15;
     577    24187341 :   const ulong CUTOFF17 = 13;
     578    24187341 :   const ulong CUTOFF18 = 11;
     579    24187341 :   const ulong CUTOFF19 = 10;
     580    24187341 :   const ulong CUTOFF20 =  9;
     581             : #else
     582     4035470 :   const ulong CUTOFF3 = 1625;
     583     4035470 :   const ulong CUTOFF4 =  255;
     584     4035470 :   const ulong CUTOFF5 =   84;
     585     4035470 :   const ulong CUTOFF6 =   40;
     586     4035470 :   const ulong CUTOFF7 =   23;
     587     4035470 :   const ulong CUTOFF8 =   15;
     588     4035470 :   const ulong CUTOFF9 =   11;
     589     4035470 :   const ulong CUTOFF10 =   9;
     590     4035470 :   const ulong CUTOFF11 =   7;
     591     4035470 :   const ulong CUTOFF12 =   6;
     592     4035470 :   const ulong CUTOFF13 =   5;
     593     4035470 :   const ulong CUTOFF14 =   4;
     594     4035470 :   const ulong CUTOFF15 =   4;
     595     4035470 :   const ulong CUTOFF16 =   3;
     596     4035470 :   const ulong CUTOFF17 =   3;
     597     4035470 :   const ulong CUTOFF18 =   3;
     598     4035470 :   const ulong CUTOFF19 =   3;
     599     4035470 :   const ulong CUTOFF20 =   3;
     600             : #endif
     601             : 
     602    28222811 :   if (p <= 2)
     603             :   {
     604     1997297 :     if (p < 2) return p;
     605     1491962 :     return k < BITS_IN_LONG? 1UL<<k: 0;
     606             :   }
     607    26225514 :   switch(k)
     608             :   {
     609             :     ulong p2, p3, p4, p5, p8;
     610     2258379 :     case 0:  return 1;
     611    10449538 :     case 1:  return p;
     612     5341660 :     case 2:  return usqru(p);
     613     2618905 :     case 3:  if (p > CUTOFF3) return 0; return p*p*p;
     614     1086148 :     case 4:  if (p > CUTOFF4) return 0; p2=p*p; return p2*p2;
     615     1238893 :     case 5:  if (p > CUTOFF5) return 0; p2=p*p; return p2*p2*p;
     616      762754 :     case 6:  if (p > CUTOFF6) return 0; p2=p*p; return p2*p2*p2;
     617      160194 :     case 7:  if (p > CUTOFF7) return 0; p2=p*p; return p2*p2*p2*p;
     618      137575 :     case 8:  if (p > CUTOFF8) return 0; p2=p*p; p4=p2*p2; return p4*p4;
     619      240380 :     case 9:  if (p > CUTOFF9) return 0; p2=p*p; p4=p2*p2; return p4*p4*p;
     620       77550 :     case 10: if (p > CUTOFF10)return 0; p2=p*p; p4=p2*p2; return p4*p4*p2;
     621       40353 :     case 11: if (p > CUTOFF11)return 0; p2=p*p; p4=p2*p2; return p4*p4*p2*p;
     622       58621 :     case 12: if (p > CUTOFF12)return 0; p2=p*p; p4=p2*p2; return p4*p4*p4;
     623       66607 :     case 13: if (p > CUTOFF13)return 0; p2=p*p; p4=p2*p2; return p4*p4*p4*p;
     624       71528 :     case 14: if (p > CUTOFF14)return 0; p2=p*p; p4=p2*p2; return p4*p4*p4*p2;
     625       89827 :     case 15: if (p > CUTOFF15)return 0;
     626       21861 :       p2=p*p; p3=p2*p; p5=p3*p2; return p5*p5*p5;
     627       63836 :     case 16: if (p > CUTOFF16)return 0;
     628       17909 :       p2=p*p; p4=p2*p2; p8=p4*p4; return p8*p8;
     629       38243 :     case 17: if (p > CUTOFF17)return 0;
     630       13245 :       p2=p*p; p4=p2*p2; p8=p4*p4; return p*p8*p8;
     631       28898 :     case 18: if (p > CUTOFF18)return 0;
     632       11495 :       p2=p*p; p4=p2*p2; p8=p4*p4; return p2*p8*p8;
     633      801961 :     case 19: if (p > CUTOFF19)return 0;
     634      788579 :       p2=p*p; p4=p2*p2; p8=p4*p4; return p*p2*p8*p8;
     635       16467 :     case 20: if (p > CUTOFF20)return 0;
     636        7346 :       p2=p*p; p4=p2*p2; p8=p4*p4; return p4*p8*p8;
     637             :   }
     638             : #ifdef LONG_IS_64BIT
     639      494199 :   switch(p)
     640             :   {
     641       56196 :     case 3: if (k > 40) return 0;
     642       24540 :       break;
     643       16206 :     case 4: if (k > 31) return 0;
     644         222 :       return 1UL<<(2*k);
     645       26961 :     case 5: if (k > 27) return 0;
     646        3915 :       break;
     647       13206 :     case 6: if (k > 24) return 0;
     648          54 :       break;
     649       18108 :     case 7: if (k > 22) return 0;
     650        1098 :       break;
     651      363522 :     default: return 0;
     652             :   }
     653             :   /* no overflow */
     654             :   {
     655       29607 :     ulong q = upowuu(p, k >> 1);
     656       29607 :     q *= q ;
     657       29607 :     return odd(k)? q*p: q;
     658             :   }
     659             : #else
     660       82998 :   return 0;
     661             : #endif
     662             : }
     663             : 
     664             : typedef struct {
     665             :   long prec, a;
     666             :   GEN (*sqr)(GEN);
     667             :   GEN (*mulug)(ulong,GEN);
     668             : } sr_muldata;
     669             : 
     670             : static GEN
     671      466693 : _rpowuu_sqr(void *data, GEN x)
     672             : {
     673      466693 :   sr_muldata *D = (sr_muldata *)data;
     674      466693 :   if (typ(x) == t_INT && lgefint(x) >= D->prec)
     675             :   { /* switch to t_REAL */
     676       17763 :     D->sqr   = &sqrr;
     677       17763 :     D->mulug = &mulur; x = itor(x, D->prec);
     678             :   }
     679      466693 :   return D->sqr(x);
     680             : }
     681             : 
     682             : static GEN
     683      171061 : _rpowuu_msqr(void *data, GEN x)
     684             : {
     685      171061 :   GEN x2 = _rpowuu_sqr(data, x);
     686      171061 :   sr_muldata *D = (sr_muldata *)data;
     687      171061 :   return D->mulug(D->a, x2);
     688             : }
     689             : 
     690             : /* return a^n as a t_REAL of precision prec. Assume a > 0, n > 0 */
     691             : GEN
     692      159105 : rpowuu(ulong a, ulong n, long prec)
     693             : {
     694             :   pari_sp av;
     695             :   GEN y, z;
     696             :   sr_muldata D;
     697             : 
     698      159105 :   if (a == 1) return real_1(prec);
     699      159105 :   if (a == 2) return real2n(n, prec);
     700      159105 :   if (n == 1) return utor(a, prec);
     701      157418 :   z = cgetr(prec);
     702      157418 :   av = avma;
     703      157418 :   D.sqr   = &sqri;
     704      157418 :   D.mulug = &mului;
     705      157418 :   D.prec = prec;
     706      157418 :   D.a = (long)a;
     707      157418 :   y = gen_powu_fold_i(utoipos(a), n, (void*)&D, &_rpowuu_sqr, &_rpowuu_msqr);
     708      157418 :   mpaff(y, z); set_avma(av); return z;
     709             : }
     710             : 
     711             : GEN
     712     7007855 : powrs(GEN x, long n)
     713             : {
     714     7007855 :   pari_sp av = avma;
     715             :   GEN y;
     716     7007855 :   if (!n) return powr0(x);
     717     7007855 :   y = gen_powu_i(x, (ulong)labs(n), NULL, &_sqrr, &_mulr);
     718     7007855 :   if (n < 0) y = invr(y);
     719     7007855 :   return gerepileuptoleaf(av,y);
     720             : }
     721             : GEN
     722      897263 : powru(GEN x, ulong n)
     723             : {
     724      897263 :   pari_sp av = avma;
     725             :   GEN y;
     726      897263 :   if (!n) return powr0(x);
     727      893959 :   y = gen_powu_i(x, n, NULL, &_sqrr, &_mulr);
     728      893959 :   return gerepileuptoleaf(av,y);
     729             : }
     730             : 
     731             : GEN
     732       14505 : powersr(GEN x, long n)
     733             : {
     734       14505 :   long prec = realprec(x);
     735       14505 :   return gen_powers(x, n, 1, &prec, &_sqrr, &_mulr, &_oner);
     736             : }
     737             : 
     738             : /* x^(s/2), assume x t_REAL */
     739             : GEN
     740           0 : powrshalf(GEN x, long s)
     741             : {
     742           0 :   if (s & 1) return sqrtr(powrs(x, s));
     743           0 :   return powrs(x, s>>1);
     744             : }
     745             : /* x^(s/2), assume x t_REAL */
     746             : GEN
     747       24487 : powruhalf(GEN x, ulong s)
     748             : {
     749       24487 :   if (s & 1) return sqrtr(powru(x, s));
     750        6335 :   return powru(x, s>>1);
     751             : }
     752             : /* x^(n/d), assume x t_REAL, return t_REAL */
     753             : GEN
     754         511 : powrfrac(GEN x, long n, long d)
     755             : {
     756             :   long z;
     757         511 :   if (!n) return powr0(x);
     758           0 :   z = cgcd(n, d); if (z > 1) { n /= z; d /= z; }
     759           0 :   if (d == 1) return powrs(x, n);
     760           0 :   x = powrs(x, n);
     761           0 :   if (d == 2) return sqrtr(x);
     762           0 :   return sqrtnr(x, d);
     763             : }
     764             : 
     765             : /* assume x != 0 */
     766             : static GEN
     767      566643 : pow_monome(GEN x, long n)
     768             : {
     769      566643 :   long i, d, dx = degpol(x);
     770             :   GEN A, b, y;
     771             : 
     772      566643 :   if (n < 0) { n = -n; y = cgetg(3, t_RFRAC); } else y = NULL;
     773             : 
     774      566643 :   if (HIGHWORD(dx) || HIGHWORD(n))
     775           8 :   {
     776             :     LOCAL_HIREMAINDER;
     777           9 :     d = (long)mulll((ulong)dx, (ulong)n);
     778           9 :     if (hiremainder || (d &~ LGBITS)) d = LGBITS; /* overflow */
     779           9 :     d += 2;
     780             :   }
     781             :   else
     782      566634 :     d = dx*n + 2;
     783      566643 :   if ((d + 1) & ~LGBITS) pari_err(e_OVERFLOW,"pow_monome [degree]");
     784      566636 :   A = cgetg(d+1, t_POL); A[1] = x[1];
     785      566636 :   for (i=2; i < d; i++) gel(A,i) = gen_0;
     786      566636 :   b = gpowgs(gel(x,dx+2), n); /* not memory clean if (n < 0) */
     787      566636 :   if (!y) y = A;
     788             :   else {
     789       20440 :     GEN c = denom_i(b);
     790       20440 :     gel(y,1) = c; if (c != gen_1) b = gmul(b,c);
     791       20440 :     gel(y,2) = A;
     792             :   }
     793      566636 :   gel(A,d) = b; return y;
     794             : }
     795             : 
     796             : /* x t_PADIC */
     797             : static GEN
     798       41517 : powps(GEN x, long n)
     799             : {
     800       41517 :   long e = n*valp(x), v;
     801       41517 :   GEN t, y, mod, p = gel(x,2);
     802             :   pari_sp av;
     803             : 
     804       41517 :   if (!signe(gel(x,4))) {
     805          84 :     if (n < 0) pari_err_INV("powps",x);
     806          77 :     return zeropadic(p, e);
     807             :   }
     808       41433 :   v = z_pval(n, p);
     809             : 
     810       41433 :   y = cgetg(5,t_PADIC);
     811       41433 :   mod = gel(x,3);
     812       41433 :   if (v == 0) mod = icopy(mod);
     813             :   else
     814             :   {
     815       40551 :     if (precp(x) == 1 && absequaliu(p, 2)) v++;
     816       40551 :     mod = mulii(mod, powiu(p,v));
     817       40551 :     mod = gerepileuptoint((pari_sp)y, mod);
     818             :   }
     819       41433 :   y[1] = evalprecp(precp(x) + v) | evalvalp(e);
     820       41433 :   gel(y,2) = icopy(p);
     821       41433 :   gel(y,3) = mod;
     822             : 
     823       41433 :   av = avma; t = gel(x,4);
     824       41433 :   if (n < 0) { t = Fp_inv(t, mod); n = -n; }
     825       41433 :   t = Fp_powu(t, n, mod);
     826       41433 :   gel(y,4) = gerepileuptoint(av, t);
     827       41433 :   return y;
     828             : }
     829             : /* x t_PADIC */
     830             : static GEN
     831         161 : powp(GEN x, GEN n)
     832             : {
     833             :   long v;
     834         161 :   GEN y, mod, p = gel(x,2);
     835             : 
     836         161 :   if (valp(x)) pari_err_OVERFLOW("valp()");
     837             : 
     838         161 :   if (!signe(gel(x,4))) {
     839          14 :     if (signe(n) < 0) pari_err_INV("powp",x);
     840           7 :     return zeropadic(p, 0);
     841             :   }
     842         147 :   v = Z_pval(n, p);
     843             : 
     844         147 :   y = cgetg(5,t_PADIC);
     845         147 :   mod = gel(x,3);
     846         147 :   if (v == 0) mod = icopy(mod);
     847             :   else
     848             :   {
     849          70 :     mod = mulii(mod, powiu(p,v));
     850          70 :     mod = gerepileuptoint((pari_sp)y, mod);
     851             :   }
     852         147 :   y[1] = evalprecp(precp(x) + v) | _evalvalp(0);
     853         147 :   gel(y,2) = icopy(p);
     854         147 :   gel(y,3) = mod;
     855         147 :   gel(y,4) = Fp_pow(gel(x,4), n, mod);
     856         147 :   return y;
     857             : }
     858             : static GEN
     859       29784 : pow_polmod(GEN x, GEN n)
     860             : {
     861       29784 :   GEN z = cgetg(3, t_POLMOD), a = gel(x,2), T = gel(x,1);
     862       29784 :   gel(z,1) = gcopy(T);
     863       29784 :   if (typ(a) != t_POL || varn(a) != varn(T) || lg(a) <= 3)
     864        3381 :     a = powgi(a, n);
     865             :   else {
     866       26403 :     pari_sp av = avma;
     867       26403 :     GEN p = NULL;
     868       26403 :     if (RgX_is_FpX(T, &p) && RgX_is_FpX(a, &p) && p)
     869             :     {
     870        7602 :       T = RgX_to_FpX(T, p); a = RgX_to_FpX(a, p);
     871        7602 :       if (lgefint(p) == 3)
     872             :       {
     873        7595 :         ulong pp = p[2];
     874        7595 :         a = Flxq_pow(ZX_to_Flx(a, pp), n, ZX_to_Flx(T, pp), pp);
     875        7595 :         a = Flx_to_ZX(a);
     876             :       }
     877             :       else
     878           7 :         a = FpXQ_pow(a, n, T, p);
     879        7602 :       a = FpX_to_mod(a, p);
     880        7602 :       a = gerepileupto(av, a);
     881             :     }
     882             :     else
     883             :     {
     884       18801 :       set_avma(av);
     885       18801 :       a = RgXQ_pow(a, n, gel(z,1));
     886             :     }
     887             :   }
     888       29784 :   gel(z,2) = a; return z;
     889             : }
     890             : 
     891             : GEN
     892   127972396 : gpowgs(GEN x, long n)
     893             : {
     894             :   long m;
     895             :   pari_sp av;
     896             :   GEN y;
     897             : 
     898   127972396 :   if (n == 0) return gpowg0(x);
     899   124427781 :   if (n == 1)
     900    72493552 :     switch (typ(x)) {
     901      802264 :       case t_QFI: return redimag(x);
     902          14 :       case t_QFR: return redreal(x);
     903    71691274 :       default: return gcopy(x);
     904             :     }
     905    51934229 :   if (n ==-1) return ginv(x);
     906    46634687 :   switch(typ(x))
     907             :   {
     908    22287113 :     case t_INT: return powis(x,n);
     909     7006847 :     case t_REAL: return powrs(x,n);
     910             :     case t_INTMOD:
     911       22239 :       y = cgetg(3,t_INTMOD); gel(y,1) = icopy(gel(x,1));
     912       22239 :       gel(y,2) = Fp_pows(gel(x,2), n, gel(x,1));
     913       22239 :       return y;
     914             :     case t_FRAC:
     915             :     {
     916      228110 :       GEN a = gel(x,1), b = gel(x,2);
     917      228110 :       long s = (signe(a) < 0 && odd(n))? -1: 1;
     918      228110 :       if (n < 0) {
     919          49 :         n = -n;
     920          49 :         if (is_pm1(a)) return powiu_sign(b, n, s); /* +-1/x[2] inverts to t_INT */
     921          42 :         swap(a, b);
     922             :       }
     923      228103 :       y = cgetg(3, t_FRAC);
     924      228103 :       gel(y,1) = powiu_sign(a, n, s);
     925      228103 :       gel(y,2) = powiu_sign(b, n, 1);
     926      228103 :       return y;
     927             :     }
     928       41517 :     case t_PADIC: return powps(x, n);
     929             :     case t_RFRAC:
     930             :     {
     931      249158 :       av = avma; y = cgetg(3, t_RFRAC); m = labs(n);
     932      249158 :       gel(y,1) = gpowgs(gel(x,1),m);
     933      249158 :       gel(y,2) = gpowgs(gel(x,2),m);
     934      249158 :       if (n < 0) y = ginv(y);
     935      249158 :       return gerepileupto(av,y);
     936             :     }
     937             :     case t_POLMOD: {
     938       29777 :       long N[] = {evaltyp(t_INT) | _evallg(3),0,0};
     939       29777 :       affsi(n,N); return pow_polmod(x, N);
     940             :     }
     941             :     case t_POL:
     942      612024 :       if (RgX_is_monomial(x)) return pow_monome(x, n);
     943             :     default: {
     944    16203283 :       pari_sp av = avma;
     945    16203283 :       y = gen_powu_i(x, (ulong)labs(n), NULL, &_sqr, &_mul);
     946    16203277 :       if (n < 0) y = ginv(y);
     947    16203277 :       return gerepileupto(av,y);
     948             :     }
     949             :   }
     950             : }
     951             : 
     952             : /* n a t_INT */
     953             : GEN
     954   114982288 : powgi(GEN x, GEN n)
     955             : {
     956             :   GEN y;
     957             : 
     958   114982288 :   if (!is_bigint(n)) return gpowgs(x, itos(n));
     959             :   /* probable overflow for non-modular types (typical exception: (X^0)^N) */
     960        8820 :   switch(typ(x))
     961             :   {
     962             :     case t_INTMOD:
     963        8527 :       y = cgetg(3,t_INTMOD); gel(y,1) = icopy(gel(x,1));
     964        8530 :       gel(y,2) = Fp_pow(gel(x,2), n, gel(x,1));
     965        8536 :       return y;
     966          59 :     case t_FFELT: return FF_pow(x,n);
     967         161 :     case t_PADIC: return powp(x, n);
     968             : 
     969             :     case t_INT:
     970          35 :       if (is_pm1(x)) return (signe(x) < 0 && mpodd(n))? gen_m1: gen_1;
     971          14 :       if (signe(x)) pari_err_OVERFLOW("lg()");
     972           7 :       if (signe(n) < 0) pari_err_INV("powgi",gen_0);
     973           7 :       return gen_0;
     974             :     case t_FRAC:
     975           7 :       pari_err_OVERFLOW("lg()");
     976             : 
     977          15 :     case t_QFR: return qfrpow(x, n);
     978           7 :     case t_POLMOD: return pow_polmod(x, n);
     979             :     default: {
     980           9 :       pari_sp av = avma;
     981           9 :       y = gen_pow_i(x, n, NULL, &_sqr, &_mul);
     982           9 :       if (signe(n) < 0) return gerepileupto(av, ginv(y));
     983           9 :       return gerepilecopy(av,y);
     984             :     }
     985             :   }
     986             : }
     987             : 
     988             : /* Assume x = 1 + O(t), n a scalar. Return x^n */
     989             : static GEN
     990        7756 : ser_pow_1(GEN x, GEN n)
     991             : {
     992             :   long lx, mi, i, j, d;
     993        7756 :   GEN y = cgetg_copy(x, &lx), X = x+2, Y = y + 2;
     994        7756 :   y[1] = evalsigne(1) | _evalvalp(0) | evalvarn(varn(x));
     995        7756 :   d = mi = lx-3; while (mi>=1 && isrationalzero(gel(X,mi))) mi--;
     996        7756 :   gel(Y,0) = gen_1;
     997      109305 :   for (i=1; i<=d; i++)
     998             :   {
     999      101549 :     pari_sp av = avma;
    1000      101549 :     GEN s = gen_0;
    1001      479178 :     for (j=1; j<=minss(i,mi); j++)
    1002             :     {
    1003      377629 :       GEN t = gsubgs(gmulgs(n,j),i-j);
    1004      377629 :       s = gadd(s, gmul(gmul(t, gel(X,j)), gel(Y,i-j)));
    1005             :     }
    1006      101549 :     gel(Y,i) = gerepileupto(av, gdivgs(s,i));
    1007             :   }
    1008        7756 :   return y;
    1009             : }
    1010             : 
    1011             : /* we suppose n != 0, valp(x) = 0 and leading-term(x) != 0. Not stack clean */
    1012             : static GEN
    1013        7875 : ser_pow(GEN x, GEN n, long prec)
    1014             : {
    1015             :   GEN y, c, lead;
    1016        7875 :   if (varncmp(gvar(n), varn(x)) <= 0) return gexp(gmul(n, glog(x,prec)), prec);
    1017        7756 :   lead = gel(x,2);
    1018        7756 :   if (gequal1(lead)) return ser_pow_1(x, n);
    1019        7434 :   x = ser_normalize(x);
    1020        7434 :   if (typ(n) == t_FRAC && !isinexact(lead) && ispower(lead, gel(n,2), &c))
    1021          77 :     c = powgi(c, gel(n,1));
    1022             :   else
    1023        7357 :     c = gpow(lead,n, prec);
    1024        7434 :   y = gmul(c, ser_pow_1(x, n));
    1025             :   /* gpow(t_POLMOD,n) can be a t_COL [conjvec] */
    1026        7434 :   if (typ(y) != t_SER) pari_err_TYPE("gpow", y);
    1027        7434 :   return y;
    1028             : }
    1029             : 
    1030             : static long
    1031        7770 : val_from_i(GEN E)
    1032             : {
    1033        7770 :   if (is_bigint(E)) pari_err_OVERFLOW("sqrtn [valuation]");
    1034        7763 :   return itos(E);
    1035             : }
    1036             : 
    1037             : /* return x^q, assume typ(x) = t_SER, typ(q) = t_INT/t_FRAC and q != 0 */
    1038             : static GEN
    1039        7777 : ser_powfrac(GEN x, GEN q, long prec)
    1040             : {
    1041        7777 :   GEN y, E = gmulsg(valp(x), q);
    1042             :   long e;
    1043             : 
    1044        7777 :   if (!signe(x))
    1045             :   {
    1046          21 :     if (gsigne(q) < 0) pari_err_INV("gpow", x);
    1047          21 :     return zeroser(varn(x), val_from_i(gfloor(E)));
    1048             :   }
    1049        7756 :   if (typ(E) != t_INT)
    1050           7 :     pari_err_DOMAIN("sqrtn", "valuation", "!=", mkintmod(gen_0, gel(q,2)), x);
    1051        7749 :   e = val_from_i(E);
    1052        7749 :   y = leafcopy(x); setvalp(y, 0);
    1053        7749 :   y = ser_pow(y, q, prec);
    1054        7749 :   setvalp(y, e); return y;
    1055             : }
    1056             : 
    1057             : static GEN
    1058         154 : gpow0(GEN x, GEN n, long prec)
    1059             : {
    1060         154 :   pari_sp av = avma;
    1061             :   long i, lx;
    1062             :   GEN y;
    1063         154 :   switch(typ(n))
    1064             :   {
    1065             :     case t_INT: case t_REAL: case t_FRAC: case t_COMPLEX: case t_QUAD:
    1066         112 :       break;
    1067             :     case t_VEC: case t_COL: case t_MAT:
    1068          35 :       y = cgetg_copy(n, &lx);
    1069          35 :       for (i=1; i<lx; i++) gel(y,i) = gpow0(x,gel(n,i),prec);
    1070          35 :       return y;
    1071           7 :     default: pari_err_TYPE("gpow(0,n)", n);
    1072             :   }
    1073         112 :   n = real_i(n);
    1074         112 :   if (gsigne(n) <= 0) pari_err_DOMAIN("gpow(0,n)", "n", "<=", gen_0, n);
    1075         105 :   if (!precision(x)) return gcopy(x);
    1076             : 
    1077          42 :   x = ground(gmulsg(gexpo(x),n));
    1078          42 :   if (is_bigint(x) || uel(x,2) >= HIGHEXPOBIT)
    1079           7 :     pari_err_OVERFLOW("gpow");
    1080          35 :   set_avma(av); return real_0_bit(itos(x));
    1081             : }
    1082             : 
    1083             : GEN
    1084    16866179 : gpow(GEN x, GEN n, long prec)
    1085             : {
    1086    16866179 :   long prec0, i, lx, tx, tn = typ(n);
    1087             :   pari_sp av;
    1088             :   GEN y;
    1089             : 
    1090    16866179 :   if (tn == t_INT) return powgi(x,n);
    1091     5039215 :   tx = typ(x);
    1092     5039215 :   if (is_matvec_t(tx))
    1093             :   {
    1094          49 :     y = cgetg_copy(x, &lx);
    1095          49 :     for (i=1; i<lx; i++) gel(y,i) = gpow(gel(x,i),n,prec);
    1096          49 :     return y;
    1097             :   }
    1098     5039166 :   av = avma;
    1099     5039166 :   switch (tx)
    1100             :   {
    1101          28 :     case t_POL: case t_RFRAC: x = toser_i(x); /* fall through */
    1102             :     case t_SER:
    1103        7574 :       if (tn == t_FRAC) return gerepileupto(av, ser_powfrac(x, n, prec));
    1104         154 :       if (valp(x))
    1105          21 :         pari_err_DOMAIN("gpow [irrational exponent]",
    1106             :                         "valuation", "!=", gen_0, x);
    1107         133 :       if (lg(x) == 2) return gerepilecopy(av, x); /* O(1) */
    1108         126 :       return gerepileupto(av, ser_pow(x, n, prec));
    1109             :   }
    1110     5031592 :   if (gequal0(x)) return gpow0(x, n, prec);
    1111     5031508 :   if (tn == t_FRAC)
    1112             :   {
    1113     4601794 :     GEN p, z, a = gel(n,1), d = gel(n,2);
    1114             :     long D;
    1115     4601794 :     switch (tx)
    1116             :     {
    1117             :     case t_INT:
    1118             :     case t_FRAC:
    1119     4654450 :       if (ispower(x, d, &z)) return powgi(z, a);
    1120       90939 :       break;
    1121             : 
    1122             :     case t_INTMOD:
    1123          21 :       p = gel(x,1);
    1124          21 :       if (!BPSW_psp(p)) pari_err_PRIME("gpow",p);
    1125          14 :       y = cgetg(3,t_INTMOD); gel(y,1) = icopy(p);
    1126          14 :       av = avma;
    1127          14 :       z = Fp_sqrtn(gel(x,2), d, p, NULL);
    1128          14 :       if (!z) pari_err_SQRTN("gpow",x);
    1129           7 :       gel(y,2) = gerepileuptoint(av, Fp_pow(z, a, p));
    1130           7 :       return y;
    1131             : 
    1132             :     case t_PADIC:
    1133          14 :       z = Qp_sqrtn(x, d, NULL); if (!z) pari_err_SQRTN("gpow",x);
    1134           7 :       return gerepileupto(av, powgi(z, a));
    1135             : 
    1136             :     case t_FFELT:
    1137          21 :       return gerepileupto(av,FF_pow(FF_sqrtn(x,d,NULL),a));
    1138             :     }
    1139     4599047 :     D = itos_or_0(d);
    1140     4599047 :     if (D == 2)
    1141             :     {
    1142     4161887 :       GEN y = gsqrt(x,prec);
    1143     4161887 :       if (!equali1(a))
    1144     4013832 :         y = gerepileupto(av, gmul(y, powgi(x, shifti(subiu(a,1), -1))));
    1145     4161887 :       return y;
    1146             :     }
    1147      437160 :     if (D && (is_real_t(tx) && gsigne(x) > 0))
    1148             :     {
    1149      396207 :       prec += nbits2extraprec(expi(a));
    1150      396207 :       if (tx != t_REAL) x = gtofp(x, prec);
    1151      396207 :       z = sqrtnr(x, D);
    1152      396207 :       if (!equali1(a)) z = powgi(z, a);
    1153      396207 :       return gerepileuptoleaf(av, z);
    1154             :     }
    1155             :   }
    1156      470667 :   i = precision(n);
    1157      470667 :   if (i) prec = i;
    1158      470667 :   prec0 = prec;
    1159      470667 :   if (!gprecision(x))
    1160             :   {
    1161       42812 :     long e = gexpo_safe(n); /* avoided if n = 0 or gexpo not defined */
    1162       42812 :     if (e > 2) prec += nbits2extraprec(e);
    1163             :   }
    1164      470667 :   y = gmul(n, glog(x,prec));
    1165      470639 :   y = gexp(y,prec);
    1166      470639 :   if (prec0 == prec) return gerepileupto(av, y);
    1167       31010 :   return gerepilecopy(av, gprec_wtrunc(y,prec0));
    1168             : }
    1169             : 
    1170             : GEN
    1171      165019 : gpowers0(GEN x, long n, GEN x0)
    1172             : {
    1173             :   long i, l;
    1174             :   GEN V;
    1175      165019 :   if (!x0) return gpowers(x,n);
    1176      152332 :   if (n < 0) return cgetg(1,t_VEC);
    1177      152332 :   l = n+2; V = cgetg(l, t_VEC); gel(V,1) = gcopy(x0);
    1178      152969 :   for (i = 2; i < l; i++) gel(V,i) = gmul(gel(V,i-1),x);
    1179      152305 :   return V;
    1180             : }
    1181             : 
    1182             : GEN
    1183      110704 : gpowers(GEN x, long n)
    1184             : {
    1185      110704 :   if (n < 0) return cgetg(1,t_VEC);
    1186      110697 :   return gen_powers(x, n, 1, (void*)x, &_sqr, &_mul, &_one);
    1187             : }
    1188             : 
    1189             : /* return [q^1,q^4,...,q^{n^2}] */
    1190             : GEN
    1191       25376 : gsqrpowers(GEN q, long n)
    1192             : {
    1193       25376 :   pari_sp av = avma;
    1194       25376 :   GEN L = gpowers0(gsqr(q), n, q); /* L[i] = q^(2i - 1), i <= n+1 */
    1195       25376 :   GEN v = cgetg(n+1, t_VEC);
    1196             :   long i;
    1197       25376 :   gel(v, 1) = gcopy(q);
    1198       25376 :   for (i = 2; i <= n ; ++i) gel(v, i) = q = gmul(q, gel(L,i)); /* q^(i^2) */
    1199       25376 :   return gerepileupto(av, v);
    1200             : }
    1201             : 
    1202             : /* 4 | N. returns a vector RU which contains exp(2*i*k*Pi/N), k=0..N-1 */
    1203             : static GEN
    1204      102885 : grootsof1_4(long N, long prec)
    1205             : {
    1206      102885 :   GEN z, RU = cgetg(N+1,t_VEC), *v  = ((GEN*)RU) + 1;
    1207      102885 :   long i, N2 = (N>>1), N4 = (N>>2), N8 = (N>>3);
    1208             :   /* z^N2 = -1, z^N4 = I; if z^k = a+I*b, then z^(N4-k) = I*conj(z) = b+a*I */
    1209             : 
    1210      102885 :   v[0] = gen_1; v[1] = z = rootsof1u_cx(N, prec);
    1211      102885 :   if (odd(N4)) N8++;
    1212      137067 :   for (i=1; i<N8; i++)
    1213             :   {
    1214       34182 :     GEN t = v[i];
    1215       34182 :     v[i+1] = gmul(z, t);
    1216       34182 :     v[N4-i] = mkcomplex(gel(t,2), gel(t,1));
    1217             :   }
    1218      102885 :   for (i=0; i<N4; i++) v[i+N4] = mulcxI(v[i]);
    1219      102885 :   for (i=0; i<N2; i++) v[i+N2] = gneg(v[i]);
    1220      102885 :   return RU;
    1221             : }
    1222             : 
    1223             : /* as above, N arbitrary */
    1224             : GEN
    1225      140591 : grootsof1(long N, long prec)
    1226             : {
    1227             :   GEN z, RU, *v;
    1228             :   long i, k;
    1229             : 
    1230      140591 :   if ((N & 3) == 0) return grootsof1_4(N, prec);
    1231       37706 :   if (N <= 2) return N == 1? mkvec(gen_1): mkvec2(gen_1, gen_m1);
    1232        9533 :   k = (N+3)>>1;
    1233        9533 :   RU = cgetg(N+1,t_VEC);
    1234        9533 :   v  = ((GEN*)RU) + 1;
    1235        9533 :   v[0] = gen_1; v[1] = z = rootsof1u_cx(N, prec);
    1236        9533 :   if (odd(N))
    1237        8066 :     for (i=2; i<k; i++) v[i] = gmul(z, v[i-1]);
    1238             :   else
    1239             :   {
    1240        1467 :     for (i=2; i<k-1; i++) v[i] = gmul(z, v[i-1]);
    1241        1467 :     v[i++] = gen_m1; /*avoid loss of accuracy*/
    1242             :   }
    1243        9533 :   for (   ; i<N; i++) v[i] = gconj(v[N-i]);
    1244        9533 :   return RU;
    1245             : }
    1246             : 
    1247             : /********************************************************************/
    1248             : /**                                                                **/
    1249             : /**                        RACINE CARREE                           **/
    1250             : /**                                                                **/
    1251             : /********************************************************************/
    1252             : /* assume x unit, e = precp(x) */
    1253             : GEN
    1254      288316 : Z2_sqrt(GEN x, long e)
    1255             : {
    1256      288316 :   ulong r = signe(x)>=0?mod16(x):16-mod16(x);
    1257             :   GEN z;
    1258             :   long ez;
    1259             :   pari_sp av;
    1260             : 
    1261      288316 :   switch(e)
    1262             :   {
    1263           7 :     case 1: return gen_1;
    1264         119 :     case 2: return (r & 3UL) == 1? gen_1: NULL;
    1265      143710 :     case 3: return (r & 7UL) == 1? gen_1: NULL;
    1266       71050 :     case 4: if (r == 1) return gen_1;
    1267       35119 :             else return (r == 9)? utoipos(3): NULL;
    1268       73430 :     default: if ((r&7UL) != 1) return NULL;
    1269             :   }
    1270       73430 :   av = avma; z = (r==1)? gen_1: utoipos(3);
    1271       73430 :   ez = 3; /* number of correct bits in z (compared to sqrt(x)) */
    1272             :   for(;;)
    1273       47978 :   {
    1274             :     GEN mod;
    1275      121408 :     ez = (ez<<1) - 1;
    1276      121408 :     if (ez > e) ez = e;
    1277      121408 :     mod = int2n(ez);
    1278      121408 :     z = addii(z, remi2n(mulii(x, Fp_inv(z,mod)), ez));
    1279      121408 :     z = shifti(z, -1); /* (z + x/z) / 2 */
    1280      121408 :     if (e == ez) return gerepileuptoint(av, z);
    1281       47978 :     if (ez < e) ez--;
    1282       47978 :     if (gc_needed(av,2))
    1283             :     {
    1284           0 :       if (DEBUGMEM > 1) pari_warn(warnmem,"Qp_sqrt");
    1285           0 :       z = gerepileuptoint(av,z);
    1286             :     }
    1287             :   }
    1288             : }
    1289             : 
    1290             : /* x unit defined modulo p^e, e > 0 */
    1291             : GEN
    1292        1764 : Qp_sqrt(GEN x)
    1293             : {
    1294        1764 :   long pp, e = valp(x);
    1295        1764 :   GEN z,y,mod, p = gel(x,2);
    1296             : 
    1297        1764 :   if (gequal0(x)) return zeropadic(p, (e+1) >> 1);
    1298        1764 :   if (e & 1) return NULL;
    1299             : 
    1300        1750 :   y = cgetg(5,t_PADIC);
    1301        1750 :   pp = precp(x);
    1302        1750 :   mod = gel(x,3);
    1303        1750 :   z   = gel(x,4); /* lift to t_INT */
    1304        1750 :   e >>= 1;
    1305        1750 :   z = Zp_sqrt(z, p, pp);
    1306        1750 :   if (!z) return NULL;
    1307        1701 :   if (absequaliu(p,2))
    1308             :   {
    1309         805 :     pp  = (pp <= 3) ? 1 : pp-1;
    1310         805 :     mod = int2n(pp);
    1311             :   }
    1312         896 :   else mod = icopy(mod);
    1313        1701 :   y[1] = evalprecp(pp) | evalvalp(e);
    1314        1701 :   gel(y,2) = icopy(p);
    1315        1701 :   gel(y,3) = mod;
    1316        1701 :   gel(y,4) = z; return y;
    1317             : }
    1318             : 
    1319             : GEN
    1320         350 : Zn_sqrt(GEN d, GEN fn)
    1321             : {
    1322         350 :   pari_sp ltop = avma, btop;
    1323         350 :   GEN b = gen_0, m = gen_1;
    1324             :   long j, np;
    1325         350 :   if (typ(d) != t_INT) pari_err_TYPE("Zn_sqrt",d);
    1326         350 :   if (typ(fn) == t_INT)
    1327           0 :     fn = absZ_factor(fn);
    1328         350 :   else if (!is_Z_factorpos(fn))
    1329           0 :     pari_err_TYPE("Zn_sqrt",fn);
    1330         350 :   np = nbrows(fn);
    1331         350 :   btop = avma;
    1332        2800 :   for (j = 1; j <= np; ++j)
    1333             :   {
    1334             :     GEN  bp, mp, pr, r;
    1335        1050 :     GEN  p = gcoeff(fn, j, 1);
    1336        1050 :     long e = itos(gcoeff(fn, j, 2));
    1337        1050 :     long v = Z_pvalrem(d,p,&r);
    1338        1050 :     if (v >= e) bp =gen_0;
    1339             :     else
    1340             :     {
    1341         952 :       if (odd(v)) return NULL;
    1342         952 :       bp = Zp_sqrt(r, p, e-v);
    1343         952 :       if (!bp)    return NULL;
    1344         952 :       if (v) bp = mulii(bp, powiu(p, v>>1L));
    1345             :     }
    1346        1050 :     mp = powiu(p, e);
    1347        1050 :     pr = mulii(m, mp);
    1348        1050 :     b = Z_chinese_coprime(b, bp, m, mp, pr);
    1349        1050 :     m = pr;
    1350        1050 :     if (gc_needed(btop, 1))
    1351           0 :       gerepileall(btop, 2, &b, &m);
    1352             :   }
    1353         350 :   return gerepileupto(ltop, b);
    1354             : }
    1355             : 
    1356             : static GEN
    1357       17094 : sqrt_ser(GEN b, long prec)
    1358             : {
    1359       17094 :   long e = valp(b), vx = varn(b), lx, lold, j;
    1360             :   ulong mask;
    1361             :   GEN a, x, lta, ltx;
    1362             : 
    1363       17094 :   if (!signe(b)) return zeroser(vx, e>>1);
    1364       17094 :   a = leafcopy(b);
    1365       17094 :   x = cgetg_copy(b, &lx);
    1366       17094 :   if (e & 1)
    1367          14 :     pari_err_DOMAIN("sqrtn", "valuation", "!=", mkintmod(gen_0, gen_2), b);
    1368       17080 :   a[1] = x[1] = evalsigne(1) | evalvarn(0) | _evalvalp(0);
    1369       17080 :   lta = gel(a,2);
    1370       17080 :   if (gequal1(lta)) ltx = lta;
    1371       14791 :   else if (!issquareall(lta,&ltx)) ltx = gsqrt(lta,prec);
    1372       17073 :   gel(x,2) = ltx;
    1373       17073 :   for (j = 3; j < lx; j++) gel(x,j) = gen_0;
    1374       17073 :   setlg(x,3);
    1375       17073 :   mask = quadratic_prec_mask(lx - 2);
    1376       17073 :   lold = 1;
    1377       99032 :   while (mask > 1)
    1378             :   {
    1379       64886 :     GEN y, x2 = gmul2n(x,1);
    1380       64886 :     long l = lold << 1, lx;
    1381             : 
    1382       64886 :     if (mask & 1) l--;
    1383       64886 :     mask >>= 1;
    1384       64886 :     setlg(a, l + 2);
    1385       64886 :     setlg(x, l + 2);
    1386       64886 :     y = sqr_ser_part(x, lold, l-1) - lold;
    1387       64886 :     for (j = lold+2; j < l+2; j++) gel(y,j) = gsub(gel(y,j), gel(a,j));
    1388       64886 :     y += lold; setvalp(y, lold);
    1389       64886 :     y = normalize(y);
    1390       64886 :     y = gsub(x, gdiv(y, x2)); /* = gmul2n(gsub(x, gdiv(a,x)), -1); */
    1391       64886 :     lx = minss(l+2, lg(y));
    1392       64886 :     for (j = lold+2; j < lx; j++) gel(x,j) = gel(y,j);
    1393       64886 :     lold = l;
    1394             :   }
    1395       17073 :   x[1] = evalsigne(1) | evalvarn(vx) | _evalvalp(e >> 1);
    1396       17073 :   return x;
    1397             : }
    1398             : 
    1399             : GEN
    1400    12776144 : gsqrt(GEN x, long prec)
    1401             : {
    1402             :   pari_sp av;
    1403             :   GEN y;
    1404             : 
    1405    12776144 :   switch(typ(x))
    1406             :   {
    1407             :     case t_INT:
    1408     3650922 :       if (!signe(x)) return real_0(prec); /* no loss of accuracy */
    1409     3650915 :       x = itor(x,prec); /* fall through */
    1410    11976649 :     case t_REAL: return sqrtr(x);
    1411             : 
    1412             :     case t_INTMOD:
    1413             :     {
    1414          35 :       GEN p = gel(x,1), a;
    1415          35 :       y = cgetg(3,t_INTMOD); gel(y,1) = icopy(p);
    1416          35 :       a = Fp_sqrt(gel(x,2),p);
    1417          21 :       if (!a)
    1418             :       {
    1419           7 :         if (!BPSW_psp(p)) pari_err_PRIME("sqrt [modulus]",p);
    1420           7 :         pari_err_SQRTN("gsqrt",x);
    1421             :       }
    1422          14 :       gel(y,2) = a; return y;
    1423             :     }
    1424             : 
    1425             :     case t_COMPLEX:
    1426             :     { /* (u+iv)^2 = a+ib <=> u^2+v^2 = sqrt(a^2+b^2), u^2-v^2=a, 2uv=b */
    1427      632153 :       GEN a = gel(x,1), b = gel(x,2), r, u, v;
    1428      632153 :       if (isrationalzero(b)) return gsqrt(a, prec);
    1429      632153 :       y = cgetg(3,t_COMPLEX); av = avma;
    1430             : 
    1431      632153 :       r = cxnorm(x);
    1432      632153 :       if (typ(r) == t_INTMOD || typ(r) == t_PADIC)
    1433           0 :         pari_err_IMPL("sqrt(complex of t_INTMODs)");
    1434      632153 :       r = gsqrt(r, prec); /* t_REAL, |a+Ib| */
    1435      632153 :       if (!signe(r))
    1436         117 :         u = v = gerepileuptoleaf(av, sqrtr(r));
    1437      632036 :       else if (gsigne(a) < 0)
    1438             :       {
    1439             :         /* v > 0 since r > 0, a < 0, rounding errors can't make the sum of two
    1440             :          * positive numbers = 0 */
    1441       66702 :         v = sqrtr( gmul2n(gsub(r,a), -1) );
    1442       66702 :         if (gsigne(b) < 0) togglesign(v);
    1443       66702 :         v = gerepileuptoleaf(av, v); av = avma;
    1444             :         /* v = 0 is impossible */
    1445       66702 :         u = gerepileuptoleaf(av, gdiv(b, shiftr(v,1)));
    1446             :       } else {
    1447      565334 :         u = sqrtr( gmul2n(gadd(r,a), -1) );
    1448      565334 :         u = gerepileuptoleaf(av, u); av = avma;
    1449      565334 :         if (!signe(u)) /* possible if a = 0.0, e.g. sqrt(0.e-10+1e-10*I) */
    1450           7 :           v = u;
    1451             :         else
    1452      565327 :           v = gerepileuptoleaf(av, gdiv(b, shiftr(u,1)));
    1453             :       }
    1454      632153 :       gel(y,1) = u;
    1455      632153 :       gel(y,2) = v; return y;
    1456             :     }
    1457             : 
    1458             :     case t_PADIC:
    1459          63 :       y = Qp_sqrt(x);
    1460          63 :       if (!y) pari_err_SQRTN("Qp_sqrt",x);
    1461          42 :       return y;
    1462             : 
    1463          70 :     case t_FFELT: return FF_sqrt(x);
    1464             : 
    1465             :     default:
    1466      167167 :       av = avma; if (!(y = toser_i(x))) break;
    1467       17094 :       return gerepilecopy(av, sqrt_ser(y, prec));
    1468             :   }
    1469      150073 :   return trans_eval("sqrt",gsqrt,x,prec);
    1470             : }
    1471             : /********************************************************************/
    1472             : /**                                                                **/
    1473             : /**                          N-th ROOT                             **/
    1474             : /**                                                                **/
    1475             : /********************************************************************/
    1476             : static void
    1477           7 : bug_logp(GEN p)
    1478             : {
    1479           7 :   if (!BPSW_psp(p)) pari_err_PRIME("p-adic log",p);
    1480           0 :   pari_err_BUG("log_p");
    1481           0 : }
    1482             : /* Let x = 1 mod p and y := (x-1)/(x+1) = 0 (p). Then
    1483             :  * log(x) = log(1+y) - log(1-y) = 2 \sum_{k odd} y^k / k.
    1484             :  * palogaux(x) returns the last sum (not multiplied by 2) */
    1485             : static GEN
    1486       25928 : palogaux(GEN x)
    1487             : {
    1488             :   long i, k, e, pp, t;
    1489       25928 :   GEN y,s,y2, p = gel(x,2);
    1490       25928 :   int is2 = absequaliu(p,2);
    1491             : 
    1492       25928 :   y = subiu(gel(x,4), 1);
    1493       25928 :   if (!signe(y))
    1494             :   {
    1495         686 :     long v = valp(x)+precp(x);
    1496         686 :     if (is2) v--;
    1497         686 :     return zeropadic(p, v);
    1498             :   }
    1499             :   /* optimize t: log(x) = log(x^(p^t)) / p^t */
    1500       25242 :   e = Z_pval(y, p); /* valp(y) = e >= 1; precp(y) = precp(x)-e */
    1501       25242 :   if (!e) bug_logp(p);
    1502       25235 :   if (is2)
    1503        2499 :     t = sqrt( (double)(precp(x)-e) / e ); /* instead of (2*e) */
    1504             :   else
    1505       22736 :     t = sqrt( (double)(precp(x)-e) / (e * (expi(p) + hammingweight(p))) );
    1506       25235 :   for (i = 0; i < t; i++) x = gpow(x, p, 0);
    1507             : 
    1508       25235 :   y = gdiv(gaddgs(x,-1), gaddgs(x,1));
    1509       25235 :   e = valp(y); /* > 0 */
    1510       25235 :   if (e <= 0) bug_logp(p);
    1511       25235 :   pp = precp(y) + e;
    1512       25235 :   if (is2) pp--;
    1513             :   else
    1514             :   {
    1515             :     GEN p1;
    1516       22736 :     for (p1=utoipos(e); abscmpui(pp,p1) > 0; pp++) p1 = mulii(p1, p);
    1517       22736 :     pp -= 2;
    1518             :   }
    1519       25235 :   k = pp/e; if (!odd(k)) k--;
    1520       25235 :   if (DEBUGLEVEL>5)
    1521           0 :     err_printf("logp: [pp,k,e,t] = [%ld,%ld,%ld,%ld]\n",pp,k,e,t);
    1522       25235 :   if (k > 1)
    1523             :   {
    1524       23513 :     y2 = gsqr(y); s = gdivgs(gen_1,k);
    1525      125349 :     while (k > 2)
    1526             :     {
    1527       78323 :       k -= 2;
    1528       78323 :       s = gadd(gmul(y2,s), gdivgs(gen_1,k));
    1529             :     }
    1530       23513 :     y = gmul(s,y);
    1531             :   }
    1532       25235 :   if (t) setvalp(y, valp(y) - t);
    1533       25235 :   return y;
    1534             : }
    1535             : 
    1536             : GEN
    1537       25935 : Qp_log(GEN x)
    1538             : {
    1539       25935 :   pari_sp av = avma;
    1540       25935 :   GEN y, p = gel(x,2), a = gel(x,4);
    1541             : 
    1542       25935 :   if (!signe(a)) pari_err_DOMAIN("Qp_log", "argument", "=", gen_0, x);
    1543       25928 :   y = leafcopy(x); setvalp(y,0);
    1544       25928 :   if (absequaliu(p,2))
    1545        2821 :     y = palogaux(gsqr(y));
    1546       23107 :   else if (gequal1(modii(a, p)))
    1547       14031 :     y = gmul2n(palogaux(y), 1);
    1548             :   else
    1549             :   { /* compute log(x^(p-1)) / (p-1) */
    1550        9076 :     GEN mod = gel(y,3), p1 = subiu(p,1);
    1551        9076 :     gel(y,4) = Fp_pow(a, p1, mod);
    1552        9076 :     p1 = diviiexact(subsi(1,mod), p1); /* 1/(p-1) */
    1553        9076 :     y = gmul(palogaux(y), shifti(p1,1));
    1554             :   }
    1555       25921 :   return gerepileupto(av,y);
    1556             : }
    1557             : 
    1558             : static GEN Qp_exp_safe(GEN x);
    1559             : 
    1560             : /*compute the p^e th root of x p-adic, assume x != 0 */
    1561             : static GEN
    1562         854 : Qp_sqrtn_ram(GEN x, long e)
    1563             : {
    1564         854 :   pari_sp ltop=avma;
    1565         854 :   GEN a, p = gel(x,2), n = powiu(p,e);
    1566         854 :   long v = valp(x), va;
    1567         854 :   if (v)
    1568             :   {
    1569             :     long z;
    1570         161 :     v = sdivsi_rem(v, n, &z);
    1571         161 :     if (z) return NULL;
    1572          91 :     x = leafcopy(x);
    1573          91 :     setvalp(x,0);
    1574             :   }
    1575             :   /*If p = 2, -1 is a root of 1 in U1: need extra check*/
    1576         784 :   if (absequaliu(p, 2) && mod8(gel(x,4)) != 1) return NULL;
    1577         749 :   a = Qp_log(x);
    1578         749 :   va = valp(a) - e;
    1579         749 :   if (va <= 0)
    1580             :   {
    1581         287 :     if (signe(gel(a,4))) return NULL;
    1582             :     /* all accuracy lost */
    1583         119 :     a = cvtop(remii(gel(x,4),p), p, 1);
    1584             :   }
    1585             :   else
    1586             :   {
    1587         462 :     setvalp(a, va); /* divide by p^e */
    1588         462 :     a = Qp_exp_safe(a);
    1589         462 :     if (!a) return NULL;
    1590             :     /* n=p^e and a^n=z*x where z is a (p-1)th-root of 1.
    1591             :      * Since z^n=z, we have (a/z)^n = x. */
    1592         462 :     a = gdiv(x, powgi(a,subiu(n,1))); /* = a/z = x/a^(n-1)*/
    1593         462 :     if (v) setvalp(a,v);
    1594             :   }
    1595         581 :   return gerepileupto(ltop,a);
    1596             : }
    1597             : 
    1598             : /*compute the nth root of x p-adic p prime with n*/
    1599             : static GEN
    1600         616 : Qp_sqrtn_unram(GEN x, GEN n, GEN *zetan)
    1601             : {
    1602             :   pari_sp av;
    1603         616 :   GEN Z, a, r, p = gel(x,2);
    1604         616 :   long v = valp(x);
    1605         616 :   if (v)
    1606             :   {
    1607             :     long z;
    1608          84 :     v = sdivsi_rem(v,n,&z);
    1609          84 :     if (z) return NULL;
    1610             :   }
    1611         609 :   r = cgetp(x); setvalp(r,v);
    1612         609 :   Z = NULL; /* -Wall */
    1613         609 :   if (zetan) Z = cgetp(x);
    1614         609 :   av = avma; a = Fp_sqrtn(gel(x,4), n, p, zetan);
    1615         609 :   if (!a) return NULL;
    1616         595 :   affii(Zp_sqrtnlift(gel(x,4), n, a, p, precp(x)), gel(r,4));
    1617         595 :   if (zetan)
    1618             :   {
    1619          14 :     affii(Zp_sqrtnlift(gen_1, n, *zetan, p, precp(x)), gel(Z,4));
    1620          14 :     *zetan = Z;
    1621             :   }
    1622         595 :   set_avma(av); return r;
    1623             : }
    1624             : 
    1625             : GEN
    1626        1183 : Qp_sqrtn(GEN x, GEN n, GEN *zetan)
    1627             : {
    1628             :   pari_sp av, tetpil;
    1629             :   GEN q, p;
    1630             :   long e;
    1631        1183 :   if (absequaliu(n, 2))
    1632             :   {
    1633          70 :     if (zetan) *zetan = gen_m1;
    1634          70 :     if (signe(n) < 0) x = ginv(x);
    1635          63 :     return Qp_sqrt(x);
    1636             :   }
    1637        1113 :   av = avma; p = gel(x,2);
    1638        1113 :   if (!signe(gel(x,4)))
    1639             :   {
    1640         203 :     if (signe(n) < 0) pari_err_INV("Qp_sqrtn", x);
    1641         203 :     q = divii(addis(n, valp(x)-1), n);
    1642         203 :     if (zetan) *zetan = gen_1;
    1643         203 :     set_avma(av); return zeropadic(p, itos(q));
    1644             :   }
    1645             :   /* treat the ramified part using logarithms */
    1646         910 :   e = Z_pvalrem(n, p, &q);
    1647         910 :   if (e) { x = Qp_sqrtn_ram(x,e); if (!x) return NULL; }
    1648         637 :   if (is_pm1(q))
    1649             :   { /* finished */
    1650          21 :     if (signe(q) < 0) x = ginv(x);
    1651          21 :     x = gerepileupto(av, x);
    1652          21 :     if (zetan)
    1653          28 :       *zetan = (e && absequaliu(p, 2))? gen_m1 /*-1 in Q_2*/
    1654          21 :                                    : gen_1;
    1655          21 :     return x;
    1656             :   }
    1657         616 :   tetpil = avma;
    1658             :   /* use hensel lift for unramified case */
    1659         616 :   x = Qp_sqrtn_unram(x, q, zetan);
    1660         616 :   if (!x) return NULL;
    1661         595 :   if (zetan)
    1662             :   {
    1663             :     GEN *gptr[2];
    1664          14 :     if (e && absequaliu(p, 2))/*-1 in Q_2*/
    1665             :     {
    1666           7 :       tetpil = avma; x = gcopy(x); *zetan = gneg(*zetan);
    1667             :     }
    1668          14 :     gptr[0] = &x; gptr[1] = zetan;
    1669          14 :     gerepilemanysp(av,tetpil,gptr,2);
    1670          14 :     return x;
    1671             :   }
    1672         581 :   return gerepile(av,tetpil,x);
    1673             : }
    1674             : 
    1675             : GEN
    1676       23927 : sqrtnint(GEN a, long n)
    1677             : {
    1678       23927 :   pari_sp ltop = avma;
    1679             :   GEN x, b, q;
    1680             :   long s, k, e;
    1681       23927 :   const ulong nm1 = n - 1;
    1682       23927 :   if (typ(a) != t_INT) pari_err_TYPE("sqrtnint",a);
    1683       23927 :   if (n <= 0) pari_err_DOMAIN("sqrtnint", "n", "<=", gen_0, stoi(n));
    1684       23920 :   if (n == 1) return icopy(a);
    1685       23920 :   s = signe(a);
    1686       23920 :   if (s < 0) pari_err_DOMAIN("sqrtnint", "x", "<", gen_0, a);
    1687       23913 :   if (!s) return gen_0;
    1688       23913 :   if (lgefint(a) == 3) return utoi(usqrtn(itou(a), n));
    1689       23504 :   e = expi(a); k = e/(2*n);
    1690       23504 :   if (k == 0)
    1691             :   {
    1692             :     long flag;
    1693         291 :     if (n > e) {set_avma(ltop); return gen_1;}
    1694         291 :     flag = cmpii(a, powuu(3, n)); set_avma(ltop);
    1695         291 :     return (flag < 0) ? gen_2: stoi(3);
    1696             :   }
    1697       23213 :   if (e < n*BITS_IN_LONG - 1)
    1698             :   {
    1699             :     ulong xs, qs;
    1700        7084 :     b = itor(a, (2*e < n*BITS_IN_LONG)? DEFAULTPREC: MEDDEFAULTPREC);
    1701        7084 :     x = mpexp(divru(logr_abs(b), n));
    1702        7084 :     xs = itou(floorr(x)) + 1; /* >= a^(1/n) */
    1703             :     for(;;) {
    1704       21216 :       q = divii(a, powuu(xs, nm1));
    1705       14150 :       if (lgefint(q) > 3) break;
    1706       14143 :       qs = itou(q); if (qs >= xs) break;
    1707        7066 :       xs -= (xs - qs + nm1)/n;
    1708             :     }
    1709        7084 :     return utoi(xs);
    1710             :   }
    1711       16129 :   b = addui(1, shifti(a, -n*k));
    1712       16129 :   x = shifti(addui(1, sqrtnint(b, n)), k);
    1713       16129 :   q = divii(a, powiu(x, nm1));
    1714       50230 :   while (cmpii(q, x) < 0) /* a priori one iteration, no GC necessary */
    1715             :   {
    1716       17972 :     x = subii(x, divis(addui(nm1, subii(x, q)), n));
    1717       17972 :     q = divii(a, powiu(x, nm1));
    1718             :   }
    1719       16129 :   return gerepileuptoleaf(ltop, x);
    1720             : }
    1721             : 
    1722             : ulong
    1723         409 : usqrtn(ulong a, ulong n)
    1724             : {
    1725             :   ulong x, s, q;
    1726         409 :   const ulong nm1 = n - 1;
    1727         409 :   if (!n) pari_err_DOMAIN("sqrtnint", "n", "=", gen_0, utoi(n));
    1728         409 :   if (n == 1 || a == 0) return a;
    1729         409 :   s = 1 + expu(a)/n; x = 1UL << s;
    1730         409 :   q = (nm1*s >= BITS_IN_LONG)? 0: a >> (nm1*s);
    1731        1622 :   while (q < x) {
    1732             :     ulong X;
    1733         804 :     x -= (x - q + nm1)/n;
    1734         804 :     X = upowuu(x, nm1);
    1735         804 :     q = X? a/X: 0;
    1736             :   }
    1737         409 :   return x;
    1738             : }
    1739             : 
    1740             : static ulong
    1741      336874 : cubic_prec_mask(long n)
    1742             : {
    1743      336874 :   long a = n, i;
    1744      336874 :   ulong mask = 0;
    1745     1927872 :   for(i = 1;; i++, mask *= 3)
    1746     1590998 :   {
    1747     1927872 :     long c = a%3;
    1748     1927872 :     if (c) mask += 3 - c;
    1749     1927872 :     a = (a+2)/3;
    1750     2264746 :     if (a==1) return mask + upowuu(3, i);
    1751             :   }
    1752             : }
    1753             : 
    1754             : /* cubic Newton iteration, |a|^(1/n), assuming a != 0 */
    1755             : GEN
    1756      465378 : sqrtnr_abs(GEN a, long n)
    1757             : {
    1758             :   pari_sp av;
    1759             :   GEN x, b;
    1760             :   long eextra, eold, n1, n2, prec, B, v;
    1761             :   ulong mask;
    1762             : 
    1763      465378 :   if (n == 1) return mpabs(a);
    1764      464819 :   if (n == 2) return sqrtr_abs(a);
    1765             : 
    1766      439412 :   prec = realprec(a);
    1767      439412 :   B = prec2nbits(prec);
    1768      439412 :   eextra = expu(n)-1;
    1769      439412 :   n1 = n+1;
    1770      439412 :   n2 = 2*n; av = avma;
    1771      439412 :   v = expo(a) / n;
    1772      439412 :   if (v) a = shiftr(a, -n*v);
    1773             : 
    1774      439412 :   b = rtor(a, DEFAULTPREC);
    1775      439412 :   x = mpexp(divru(logr_abs(b), n));
    1776      439412 :   if (prec == DEFAULTPREC)
    1777             :   {
    1778      121149 :     if (v) shiftr_inplace(x, v);
    1779      121149 :     return gerepileuptoleaf(av, x);
    1780             :   }
    1781      318263 :   mask = cubic_prec_mask(B + 63);
    1782      318263 :   eold = 1;
    1783             :   for(;;)
    1784     1257421 :   { /* reach 64 */
    1785     1575684 :     long enew = eold * 3;
    1786     1575684 :     enew -= mask % 3;
    1787     1575684 :     if (enew > 64) break; /* back up one step */
    1788     1257421 :     mask /= 3;
    1789     1257421 :     eold = enew;
    1790             :   }
    1791             :   for(;;)
    1792      248770 :   {
    1793      567033 :     long pr, enew = eold * 3;
    1794             :     GEN y, z;
    1795      567033 :     enew -= mask % 3;
    1796      567033 :     mask /= 3;
    1797      567033 :     pr = nbits2prec(enew + eextra);
    1798      567033 :     b = rtor(a, pr); setsigne(b,1);
    1799      567033 :     x = rtor(x, pr);
    1800      567033 :     y = subrr(powru(x, n), b);
    1801      567033 :     z = divrr(y, addrr(mulur(n1, y), mulur(n2, b)));
    1802      567033 :     shiftr_inplace(z,1);
    1803      567033 :     x = mulrr(x, subsr(1,z));
    1804      567033 :     if (mask == 1)
    1805             :     {
    1806      318263 :       if (v) shiftr_inplace(x, v);
    1807      318263 :       return gerepileuptoleaf(av, gprec_wtrunc(x,prec));
    1808             :     }
    1809      248770 :     eold = enew;
    1810             :   }
    1811             : }
    1812             : 
    1813             : static void
    1814       31036 : shiftc_inplace(GEN z, long d)
    1815             : {
    1816       31036 :   shiftr_inplace(gel(z,1), d);
    1817       31036 :   shiftr_inplace(gel(z,2), d);
    1818       31036 : }
    1819             : 
    1820             : /* exp(2*Pi*I/n), same iteration as sqrtnr_abs, different initial point */
    1821             : static GEN
    1822      103833 : sqrtnof1(ulong n, long prec)
    1823             : {
    1824             :   pari_sp av;
    1825             :   GEN x;
    1826             :   long eold, n1, n2, B;
    1827             :   ulong mask;
    1828             : 
    1829      103833 :   B = prec2nbits(prec);
    1830      103833 :   n1 = n+1;
    1831      103833 :   n2 = 2*n; av = avma;
    1832             : 
    1833      103833 :   x = expIr(divru(Pi2n(1, LOWDEFAULTPREC), n));
    1834      103833 :   if (prec == LOWDEFAULTPREC) return gerepileupto(av, x);
    1835       18611 :   mask = cubic_prec_mask(B + BITS_IN_LONG-1);
    1836       18611 :   eold = 1;
    1837             :   for(;;)
    1838       72382 :   { /* reach BITS_IN_LONG */
    1839       90993 :     long enew = eold * 3;
    1840       90993 :     enew -= mask % 3;
    1841       90993 :     if (enew > BITS_IN_LONG) break; /* back up one step */
    1842       72382 :     mask /= 3;
    1843       72382 :     eold = enew;
    1844             :   }
    1845             :   for(;;)
    1846       12425 :   {
    1847       31036 :     long pr, enew = eold * 3;
    1848             :     GEN y, z;
    1849       31036 :     enew -= mask % 3;
    1850       31036 :     mask /= 3;
    1851       31036 :     pr = nbits2prec(enew);
    1852       31036 :     x = cxtofp(x, pr);
    1853       31036 :     y = gsub(gpowgs(x, n), gen_1);
    1854       31036 :     z = gdiv(y, gaddgs(gmulsg(n1, y), n2));
    1855       31036 :     shiftc_inplace(z,1);
    1856       31036 :     x = gmul(x, gsubsg(1, z));
    1857       31036 :     if (mask == 1) return gerepilecopy(av, gprec_w(x,prec));
    1858       12425 :     eold = enew;
    1859             :   }
    1860             : }
    1861             : 
    1862             : /* exp(2iPi/d) */
    1863             : GEN
    1864      238829 : rootsof1u_cx(ulong n, long prec)
    1865             : {
    1866      238829 :   switch(n)
    1867             :   {
    1868       11718 :     case 1: return gen_1;
    1869        2548 :     case 2: return gen_m1;
    1870       41845 :     case 4: return gen_I();
    1871             :     case 3: case 6: case 12:
    1872             :     {
    1873        6111 :       pari_sp av = avma;
    1874        6111 :       GEN a = (n == 3)? mkfrac(gen_m1,gen_2): ghalf;
    1875        6111 :       GEN sq3 = sqrtr_abs(utor(3, prec));
    1876        6111 :       shiftr_inplace(sq3, -1);
    1877        6111 :       a = (n == 12)? mkcomplex(sq3, a): mkcomplex(a, sq3);
    1878        6111 :       return gerepilecopy(av, a);
    1879             :     }
    1880             :     case 8:
    1881             :     {
    1882       72774 :       pari_sp av = avma;
    1883       72774 :       GEN sq2 = sqrtr_abs(utor(2, prec));
    1884       72774 :       shiftr_inplace(sq2,-1);
    1885       72774 :       return gerepilecopy(av, mkcomplex(sq2, sq2));
    1886             :     }
    1887             :   }
    1888      103833 :   return sqrtnof1(n, prec);
    1889             : }
    1890             : /* e(a/b) */
    1891             : GEN
    1892       12782 : rootsof1q_cx(long a, long b, long prec)
    1893             : {
    1894       12782 :   long g = cgcd(a,b);
    1895             :   GEN z;
    1896       12782 :   if (g != 1) { a /= g; b /= g; }
    1897       12782 :   if (b < 0) { b = -b; a = -a; }
    1898       12782 :   z = rootsof1u_cx(b, prec);
    1899       12782 :   if (a < 0) { z = conj_i(z); a = -a; }
    1900       12782 :   return gpowgs(z, a);
    1901             : }
    1902             : 
    1903             : /* initializes powers of e(a/b) */
    1904             : GEN
    1905       12964 : rootsof1powinit(long a, long b, long prec)
    1906             : {
    1907       12964 :   long g = cgcd(a,b);
    1908       12964 :   if (g != 1) { a /= g; b /= g; }
    1909       12964 :   if (b < 0) { b = -b; a = -a; }
    1910       12964 :   a %= b; if (a < 0) a += b;
    1911       12964 :   return mkvec2(grootsof1(b,prec), mkvecsmall2(a,b));
    1912             : }
    1913             : /* T = rootsof1powinit(a,b); return  e(a/b)^c */
    1914             : GEN
    1915    12301716 : rootsof1pow(GEN T, long c)
    1916             : {
    1917    12301716 :   GEN vz = gel(T,1), ab = gel(T,2);
    1918    12301716 :   long a = ab[1], b = ab[2]; /* a >= 0, b > 0 */
    1919    12301716 :   c %= b; if (c < 0) c += b;
    1920    12301716 :   a = Fl_mul(a, c, b);
    1921    12301716 :   return gel(vz, a + 1);
    1922             : }
    1923             : 
    1924             : /* exp(2iPi/d), assume d a t_INT */
    1925             : GEN
    1926        3780 : rootsof1_cx(GEN d, long prec)
    1927             : {
    1928        3780 :   if (lgefint(d) == 3) return rootsof1u_cx((ulong)d[2], prec);
    1929           0 :   return expIr(divri(Pi2n(1,prec), d));
    1930             : }
    1931             : 
    1932             : GEN
    1933        3774 : gsqrtn(GEN x, GEN n, GEN *zetan, long prec)
    1934             : {
    1935             :   long i, lx, tx;
    1936             :   pari_sp av;
    1937             :   GEN y, z;
    1938        3774 :   if (typ(n)!=t_INT) pari_err_TYPE("sqrtn",n);
    1939        3774 :   if (!signe(n)) pari_err_DOMAIN("sqrtn", "n", "=", gen_0, n);
    1940        3774 :   if (is_pm1(n))
    1941             :   {
    1942          70 :     if (zetan) *zetan = gen_1;
    1943          70 :     return (signe(n) > 0)? gcopy(x): ginv(x);
    1944             :   }
    1945        3704 :   if (zetan) *zetan = gen_0;
    1946        3704 :   tx = typ(x);
    1947        3704 :   if (is_matvec_t(tx))
    1948             :   {
    1949           7 :     y = cgetg_copy(x, &lx);
    1950           7 :     for (i=1; i<lx; i++) gel(y,i) = gsqrtn(gel(x,i),n,NULL,prec);
    1951           7 :     return y;
    1952             :   }
    1953        3697 :   av = avma;
    1954        3697 :   switch(tx)
    1955             :   {
    1956             :   case t_INTMOD:
    1957             :     {
    1958         175 :       GEN p = gel(x,1), s;
    1959         175 :       z = gen_0;
    1960         175 :       y = cgetg(3,t_INTMOD);  gel(y,1) = icopy(p);
    1961         175 :       if (zetan) { z = cgetg(3,t_INTMOD); gel(z,1) = gel(y,1); }
    1962         175 :       s = Fp_sqrtn(gel(x,2),n,p,zetan);
    1963         154 :       if (!s) {
    1964          28 :         if (zetan) {set_avma(av); return gen_0;}
    1965          21 :         if (!BPSW_psp(p)) pari_err_PRIME("sqrtn [modulus]",p);
    1966          14 :         pari_err_SQRTN("gsqrtn",x);
    1967             :       }
    1968         126 :       gel(y,2) = s;
    1969         126 :       if (zetan) { gel(z,2) = *zetan; *zetan = z; }
    1970         126 :       return y;
    1971             :     }
    1972             : 
    1973             :   case t_PADIC:
    1974          56 :     y = Qp_sqrtn(x,n,zetan);
    1975          49 :     if (!y) {
    1976           7 :       if (zetan) return gen_0;
    1977           7 :       pari_err_SQRTN("gsqrtn",x);
    1978             :     }
    1979          42 :     return y;
    1980             : 
    1981         196 :   case t_FFELT: return FF_sqrtn(x,n,zetan);
    1982             : 
    1983             :   case t_INT: case t_FRAC: case t_REAL: case t_COMPLEX:
    1984        2906 :     i = precision(x); if (i) prec = i;
    1985        2906 :     if (isint1(x))
    1986           7 :       y = real_1(prec);
    1987        2899 :     else if (gequal0(x))
    1988             :     {
    1989             :       long b;
    1990          21 :       if (signe(n) < 0) pari_err_INV("gsqrtn",x);
    1991          21 :       if (isinexactreal(x))
    1992          14 :         b = sdivsi(gexpo(x), n);
    1993             :       else
    1994           7 :         b = -prec2nbits(prec);
    1995          21 :       if (typ(x) == t_COMPLEX)
    1996             :       {
    1997           7 :         y = cgetg(3,t_COMPLEX);
    1998           7 :         gel(y,1) = gel(y,2) = real_0_bit(b);
    1999             :       }
    2000             :       else
    2001          14 :         y = real_0_bit(b);
    2002             :     }
    2003             :     else
    2004             :     {
    2005        2878 :       long nn = itos_or_0(n);
    2006        2878 :       if (tx == t_INT) { x = itor(x,prec); tx = t_REAL; }
    2007        2878 :       if (nn > 0 && tx == t_REAL && signe(x) > 0)
    2008        1103 :         y = sqrtnr(x, nn);
    2009             :       else
    2010        1775 :         y = gexp(gdiv(glog(x,prec), n), prec);
    2011        2878 :       y = gerepileupto(av, y);
    2012             :     }
    2013        2906 :     if (zetan) *zetan = rootsof1_cx(n, prec);
    2014        2906 :     return y;
    2015             : 
    2016             :   case t_QUAD:
    2017           7 :     return gsqrtn(quadtofp(x, prec), n, zetan, prec);
    2018             : 
    2019             :   default:
    2020         357 :     av = avma; if (!(y = toser_i(x))) break;
    2021         357 :     return gerepileupto(av, ser_powfrac(y, ginv(n), prec));
    2022             :   }
    2023           0 :   pari_err_TYPE("sqrtn",x);
    2024             :   return NULL;/* LCOV_EXCL_LINE */
    2025             : }
    2026             : 
    2027             : /********************************************************************/
    2028             : /**                                                                **/
    2029             : /**                             EXP(X) - 1                         **/
    2030             : /**                                                                **/
    2031             : /********************************************************************/
    2032             : /* exp(|x|) - 1, assume x != 0.
    2033             :  * For efficiency, x should be reduced mod log(2): if so, we have a < 0 */
    2034             : GEN
    2035    10391130 : exp1r_abs(GEN x)
    2036             : {
    2037    10391130 :   long l = realprec(x), a = expo(x), b = prec2nbits(l), L, i, n, m, B;
    2038             :   GEN y, p2, X;
    2039             :   pari_sp av;
    2040             :   double d;
    2041             : 
    2042    10391130 :   if (b + a <= 0) return mpabs(x);
    2043             : 
    2044    10375492 :   y = cgetr(l); av = avma;
    2045    10375492 :   B = b/3 + BITS_IN_LONG + (BITS_IN_LONG*BITS_IN_LONG)/ b;
    2046    10375492 :   d = a/2.; m = (long)(d + sqrt(d*d + B)); /* >= 0 */
    2047    10375492 :   if (m < (-a) * 0.1) m = 0; /* not worth it */
    2048    10375492 :   L = l + nbits2extraprec(m);
    2049             :  /* Multiplication is quadratic in this range (l is small, otherwise we
    2050             :   * use logAGM + Newton). Set Y = 2^(-e-a) x, compute truncated series
    2051             :   * sum x^k/k!: this costs roughly
    2052             :   *    m b^2 + sum_{k <= n} (k e + BITS_IN_LONG)^2
    2053             :   * bit operations with |x| <  2^(1+a), |Y| < 2^(1-e) , m = e+a and b bits of
    2054             :   * accuracy needed, so
    2055             :   *    B := (b / 3 + BITS_IN_LONG + BITS_IN_LONG^2 / b) ~ m(m-a)
    2056             :   * we want b ~ 3 m (m-a) or m~b+a hence
    2057             :   *     m = min( a/2 + sqrt(a^2/4 + B),  b + a )
    2058             :   * NB: e ~ (b/3)^(1/2) as b -> oo
    2059             :   *
    2060             :   * Truncate the sum at k = n (>= 1), the remainder is
    2061             :   *   sum_{k >= n+1} Y^k / k! < Y^(n+1) / (n+1)! (1-Y) < Y^(n+1) / n!
    2062             :   * We want Y^(n+1) / n! <= Y 2^-b, hence -n log_2 |Y| + log_2 n! >= b
    2063             :   *   log n! ~ (n + 1/2) log(n+1) - (n+1) + log(2Pi)/2,
    2064             :   * error bounded by 1/6(n+1) <= 1/12. Finally, we want
    2065             :   * n (-1/log(2) -log_2 |Y| + log_2(n+1)) >= b  */
    2066    10375492 :   b += m;
    2067    10375492 :   d = m-dbllog2(x)-1/M_LN2; /* ~ -log_2 Y - 1/log(2) */
    2068    10375492 :   n = (long)(b / d);
    2069    10375492 :   if (n > 1)
    2070    10359026 :     n = (long)(b / (d + log2((double)n+1))); /* log~constant in small ranges */
    2071    10375492 :   while (n*(d+log2((double)n+1)) < b) n++; /* expect few corrections */
    2072             : 
    2073    10375492 :   X = rtor(x,L); shiftr_inplace(X, -m); setsigne(X, 1);
    2074    10375492 :   if (n == 1) p2 = X;
    2075             :   else
    2076             :   {
    2077    10375492 :     long s = 0, l1 = nbits2prec((long)(d + n + 16));
    2078    10375492 :     GEN unr = real_1(L);
    2079             :     pari_sp av2;
    2080             : 
    2081    10375492 :     p2 = cgetr(L); av2 = avma;
    2082   146055644 :     for (i=n; i>=2; i--, set_avma(av2))
    2083             :     { /* compute X^(n-1)/n! + ... + X/2 + 1 */
    2084             :       GEN p1, p3;
    2085   135680152 :       setprec(X,l1); p3 = divru(X,i);
    2086   135680152 :       l1 += dvmdsBIL(s - expo(p3), &s); if (l1>L) l1=L;
    2087   135680152 :       setprec(unr,l1); p1 = addrr_sign(unr,1, i == n? p3: mulrr(p3,p2),1);
    2088   135680152 :       setprec(p2,l1); affrr(p1,p2); /* p2 <- 1 + (X/i)*p2 */
    2089             :     }
    2090    10375492 :     setprec(X,L); p2 = mulrr(X,p2);
    2091             :   }
    2092             : 
    2093   113787178 :   for (i=1; i<=m; i++)
    2094             :   {
    2095   103411686 :     if (realprec(p2) > L) setprec(p2,L);
    2096   103411686 :     p2 = mulrr(p2, addsr(2,p2));
    2097             :   }
    2098    10375492 :   affrr_fixlg(p2,y); set_avma(av); return y;
    2099             : }
    2100             : 
    2101             : GEN
    2102        9026 : mpexpm1(GEN x)
    2103             : {
    2104        9026 :   const long s = 6;
    2105        9026 :   long l, sx = signe(x);
    2106             :   GEN y, z;
    2107             :   pari_sp av;
    2108        9026 :   if (!sx) return real_0_bit(expo(x));
    2109        9019 :   l = realprec(x);
    2110        9019 :   if (l > maxss(EXPNEWTON_LIMIT, (1L<<s) + 2))
    2111             :   {
    2112           6 :     long e = expo(x);
    2113           6 :     if (e < 0) x = rtor(x, l + nbits2extraprec(-e));
    2114           6 :     return subrs(mpexp(x), 1);
    2115             :   }
    2116        9013 :   if (sx > 0) return exp1r_abs(x);
    2117             :   /* compute exp(x) * (1 - exp(-x)) */
    2118        4156 :   av = avma; y = exp1r_abs(x);
    2119        4156 :   z = addsr(1, y); setsigne(z, -1);
    2120        4156 :   return gerepileupto(av, divrr(y, z));
    2121             : }
    2122             : 
    2123             : static GEN serexp(GEN x, long prec);
    2124             : GEN
    2125       10615 : gexpm1(GEN x, long prec)
    2126             : {
    2127       10615 :   switch(typ(x))
    2128             :   {
    2129        4088 :     case t_REAL: return mpexpm1(x);
    2130        4994 :     case t_COMPLEX: return cxexpm1(x,prec);
    2131          14 :     case t_PADIC: return gsubgs(Qp_exp(x), 1);
    2132             :     default:
    2133             :     {
    2134        1519 :       pari_sp av = avma;
    2135             :       long ey;
    2136             :       GEN y;
    2137        1519 :       if (!(y = toser_i(x))) break;
    2138        1498 :       ey = valp(y);
    2139        1498 :       if (ey < 0) pari_err_DOMAIN("expm1","valuation", "<", gen_0, x);
    2140        1498 :       if (gequal0(y)) return gcopy(y);
    2141        1491 :       if (ey)
    2142         343 :         return gerepileupto(av, gsubgs(serexp(y,prec), 1));
    2143             :       else
    2144             :       {
    2145        1148 :         GEN e1 = gexpm1(gel(y,2), prec), e = gaddgs(e1,1);
    2146        1148 :         y = gmul(e, serexp(serchop0(y),prec));
    2147        1148 :         gel(y,2) = e1;
    2148        1148 :         return gerepilecopy(av, y);
    2149             :       }
    2150             :     }
    2151             :   }
    2152          21 :   return trans_eval("expm1",gexpm1,x,prec);
    2153             : }
    2154             : /********************************************************************/
    2155             : /**                                                                **/
    2156             : /**                             EXP(X)                             **/
    2157             : /**                                                                **/
    2158             : /********************************************************************/
    2159             : 
    2160             : /* centermod(x, log(2)), set *sh to the quotient */
    2161             : static GEN
    2162    10330605 : modlog2(GEN x, long *sh)
    2163             : {
    2164    10330605 :   double d = rtodbl(x), da = fabs(d);
    2165    10330605 :   long q = (long) ((da + (M_LN2/2))/M_LN2);
    2166    10330605 :   if (da > M_LN2 * LONG_MAX)
    2167          14 :     pari_err_OVERFLOW("expo()"); /* avoid overflow in  q */
    2168    10330591 :   if (d < 0) q = -q;
    2169    10330591 :   *sh = q;
    2170    10330591 :   if (q) {
    2171     9810057 :     long l = realprec(x) + 1;
    2172     9810057 :     x = subrr(rtor(x,l), mulsr(q, mplog2(l)));
    2173     9810057 :     if (!signe(x)) return NULL;
    2174             :   }
    2175    10330591 :   return x;
    2176             : }
    2177             : 
    2178             : static GEN
    2179    10329908 : mpexp_basecase(GEN x)
    2180             : {
    2181    10329908 :   pari_sp av = avma;
    2182    10329908 :   long sh, l = realprec(x);
    2183             :   GEN y, z;
    2184             : 
    2185    10329908 :   y = modlog2(x, &sh);
    2186    10329894 :   if (!y) { set_avma(av); return real2n(sh, l); }
    2187    10329894 :   z = addsr(1, exp1r_abs(y));
    2188    10329894 :   if (signe(y) < 0) z = invr(z);
    2189    10329894 :   if (sh) {
    2190     9809366 :     shiftr_inplace(z, sh);
    2191     9809366 :     if (realprec(z) > l) z = rtor(z, l); /* spurious precision increase */
    2192             :   }
    2193             : #ifdef DEBUG
    2194             : {
    2195             :   GEN t = mplog(z), u = divrr(subrr(x, t),x);
    2196             :   if (signe(u) && expo(u) > 5-prec2nbits(minss(l,realprec(t))))
    2197             :     pari_err_BUG("exp");
    2198             : }
    2199             : #endif
    2200    10329894 :   return gerepileuptoleaf(av, z); /* NOT affrr, precision often increases */
    2201             : }
    2202             : 
    2203             : GEN
    2204    10542307 : mpexp(GEN x)
    2205             : {
    2206    10542307 :   const long s = 6; /*Initial steps using basecase*/
    2207    10542307 :   long i, p, l = realprec(x), sh;
    2208             :   GEN a, t, z;
    2209             :   ulong mask;
    2210             : 
    2211    10542307 :   if (l <= maxss(EXPNEWTON_LIMIT, (1L<<s) + 2))
    2212             :   {
    2213    10541610 :     if (!signe(x)) return mpexp0(x);
    2214    10329211 :     return mpexp_basecase(x);
    2215             :   }
    2216         697 :   z = cgetr(l); /* room for result */
    2217         697 :   x = modlog2(x, &sh);
    2218         697 :   if (!x) { set_avma((pari_sp)(z+lg(z))); return real2n(sh, l); }
    2219         697 :   constpi(l); /* precompute for later logr_abs() */
    2220         697 :   mask = quadratic_prec_mask(prec2nbits(l)+BITS_IN_LONG);
    2221         697 :   for(i=0, p=1; i<s+TWOPOTBITS_IN_LONG; i++) { p <<= 1; if (mask & 1) p-=1; mask >>= 1; }
    2222         697 :   a = mpexp_basecase(rtor(x, nbits2prec(p)));
    2223         697 :   x = addrs(x,1);
    2224         697 :   if (realprec(x) < l+EXTRAPRECWORD) x = rtor(x, l+EXTRAPRECWORD);
    2225         697 :   a = rtor(a, l+EXTRAPRECWORD); /*append 0s */
    2226         697 :   t = NULL;
    2227             :   for(;;)
    2228             :   {
    2229         771 :     p <<= 1; if (mask & 1) p--;
    2230         734 :     mask >>= 1;
    2231         734 :     setprec(x, nbits2prec(p));
    2232         734 :     setprec(a, nbits2prec(p));
    2233         734 :     t = mulrr(a, subrr(x, logr_abs(a))); /* a (x - log(a)) */
    2234         734 :     if (mask == 1) break;
    2235          37 :     affrr(t, a); set_avma((pari_sp)a);
    2236             :   }
    2237         697 :   affrr(t,z);
    2238         697 :   if (sh) shiftr_inplace(z, sh);
    2239         697 :   set_avma((pari_sp)z); return z;
    2240             : }
    2241             : 
    2242             : static long
    2243       20146 : Qp_exp_prec(GEN x)
    2244             : {
    2245       20146 :   long k, e = valp(x), n = e + precp(x);
    2246       20146 :   GEN p = gel(x,2);
    2247       20146 :   int is2 = absequaliu(p,2);
    2248       20146 :   if (e < 1 || (e == 1 && is2)) return -1;
    2249       20118 :   if (is2)
    2250             :   {
    2251        6321 :     n--; e--; k = n/e;
    2252        6321 :     if (n%e == 0) k--;
    2253             :   }
    2254             :   else
    2255             :   { /* e > 0, n > 0 */
    2256       13797 :     GEN r, t = subiu(p, 1);
    2257       13797 :     k = itos(dvmdii(subiu(muliu(t,n), 1), subiu(muliu(t,e), 1), &r));
    2258       13797 :     if (!signe(r)) k--;
    2259             :   }
    2260       20118 :   return k;
    2261             : }
    2262             : 
    2263             : static GEN
    2264       21574 : Qp_exp_safe(GEN x)
    2265             : {
    2266             :   long k;
    2267             :   pari_sp av;
    2268             :   GEN y;
    2269             : 
    2270       21574 :   if (gequal0(x)) return gaddgs(x,1);
    2271       20076 :   k = Qp_exp_prec(x);
    2272       20076 :   if (k < 0) return NULL;
    2273       20069 :   av = avma;
    2274       20069 :   for (y=gen_1; k; k--) y = gaddsg(1, gdivgs(gmul(y,x), k));
    2275       20069 :   return gerepileupto(av, y);
    2276             : }
    2277             : 
    2278             : GEN
    2279       21112 : Qp_exp(GEN x)
    2280             : {
    2281       21112 :   GEN y = Qp_exp_safe(x);
    2282       21112 :   if (!y) pari_err_DOMAIN("gexp(t_PADIC)","argument","",gen_0,x);
    2283       21105 :   return y;
    2284             : }
    2285             : 
    2286             : static GEN
    2287          49 : cos_p(GEN x)
    2288             : {
    2289             :   long k;
    2290             :   pari_sp av;
    2291             :   GEN x2, y;
    2292             : 
    2293          49 :   if (gequal0(x)) return gaddgs(x,1);
    2294          28 :   k = Qp_exp_prec(x);
    2295          28 :   if (k < 0) return NULL;
    2296          21 :   av = avma; x2 = gsqr(x);
    2297          21 :   if (k & 1) k--;
    2298         105 :   for (y=gen_1; k; k-=2)
    2299             :   {
    2300          84 :     GEN t = gdiv(gmul(y,x2), muluu(k, k-1));
    2301          84 :     y = gsubsg(1, t);
    2302             :   }
    2303          21 :   return gerepileupto(av, y);
    2304             : }
    2305             : static GEN
    2306          63 : sin_p(GEN x)
    2307             : {
    2308             :   long k;
    2309             :   pari_sp av;
    2310             :   GEN x2, y;
    2311             : 
    2312          63 :   if (gequal0(x)) return gcopy(x);
    2313          42 :   k = Qp_exp_prec(x);
    2314          42 :   if (k < 0) return NULL;
    2315          28 :   av = avma; x2 = gsqr(x);
    2316          28 :   if (k & 1) k--;
    2317         133 :   for (y=gen_1; k; k-=2)
    2318             :   {
    2319         105 :     GEN t = gdiv(gmul(y,x2), muluu(k, k+1));
    2320         105 :     y = gsubsg(1, t);
    2321             :   }
    2322          28 :   return gerepileupto(av, gmul(y, x));
    2323             : }
    2324             : 
    2325             : static GEN
    2326     1414313 : cxexp(GEN x, long prec)
    2327             : {
    2328     1414313 :   GEN r, p1, p2, y = cgetg(3,t_COMPLEX);
    2329     1414313 :   pari_sp av = avma, tetpil;
    2330             :   long l;
    2331     1414313 :   l = precision(x); if (l > prec) prec = l;
    2332     1414313 :   r = gexp(gel(x,1),prec);
    2333     1414313 :   if (gequal0(r)) { gel(y,1) = r; gel(y,2) = r; return y; }
    2334     1414313 :   gsincos(gel(x,2),&p2,&p1,prec);
    2335     1414313 :   tetpil = avma;
    2336     1414313 :   gel(y,1) = gmul(r,p1);
    2337     1414313 :   gel(y,2) = gmul(r,p2);
    2338     1414313 :   gerepilecoeffssp(av,tetpil,y+1,2);
    2339     1414313 :   return y;
    2340             : }
    2341             : 
    2342             : /* given a t_SER x^v s(x), with s(0) != 0, return x^v(s - s(0)), shallow */
    2343             : GEN
    2344       28973 : serchop0(GEN s)
    2345             : {
    2346       28973 :   long i, l = lg(s);
    2347             :   GEN y;
    2348       28973 :   if (l == 2) return s;
    2349       28973 :   if (l == 3 && isexactzero(gel(s,2))) return s;
    2350       28973 :   y = cgetg(l, t_SER); y[1] = s[1];
    2351       28973 :   gel(y,2) = gen_0; for (i=3; i <l; i++) gel(y,i) = gel(s,i);
    2352       28973 :   return normalize(y);
    2353             : }
    2354             : 
    2355             : GEN
    2356          42 : serchop_i(GEN s, long n)
    2357             : {
    2358          42 :   long i, m, l = lg(s);
    2359             :   GEN y;
    2360          42 :   if (l == 2 || (l == 3 && isexactzero(gel(s,2))))
    2361             :   {
    2362          14 :     if (valp(s) < n) { s = shallowcopy(s); setvalp(s,n); }
    2363          14 :     return s;
    2364             :   }
    2365          28 :   m = n - valp(s); if (m < 0) return s;
    2366          21 :   if (l-m <= 2) return zeroser(varn(s), n);
    2367          14 :   y = cgetg(l-m, t_SER); y[1] = s[1]; setvalp(y, valp(y)+m);
    2368          14 :   for (i=m+2; i < l; i++) gel(y,i-m) = gel(s,i);
    2369          14 :   return normalize(y);
    2370             : }
    2371             : GEN
    2372          42 : serchop(GEN s, long n)
    2373             : {
    2374          42 :   pari_sp av = avma;
    2375          42 :   if (typ(s) != t_SER) pari_err_TYPE("serchop",s);
    2376          42 :   return gerepilecopy(av, serchop_i(s,n));
    2377             : }
    2378             : 
    2379             : static GEN
    2380       76419 : serexp(GEN x, long prec)
    2381             : {
    2382             :   pari_sp av;
    2383             :   long i,j,lx,ly,ex,mi;
    2384             :   GEN p1,y,xd,yd;
    2385             : 
    2386       76419 :   ex = valp(x);
    2387       76419 :   if (ex < 0) pari_err_DOMAIN("exp","valuation", "<", gen_0, x);
    2388       76412 :   if (gequal0(x)) return gaddsg(1,x);
    2389       56077 :   lx = lg(x);
    2390       56077 :   if (ex)
    2391             :   {
    2392       40005 :     ly = lx+ex; y = cgetg(ly,t_SER);
    2393       40005 :     mi = lx-1; while (mi>=3 && isrationalzero(gel(x,mi))) mi--;
    2394       40005 :     mi += ex-2;
    2395       40005 :     y[1] = evalsigne(1) | _evalvalp(0) | evalvarn(varn(x));
    2396             :     /* zd[i] = coefficient of X^i in z */
    2397       40005 :     xd = x+2-ex; yd = y+2; ly -= 2;
    2398       40005 :     gel(yd,0) = gen_1;
    2399       40005 :     for (i=1; i<ex; i++) gel(yd,i) = gen_0;
    2400      422856 :     for (   ; i<ly; i++)
    2401             :     {
    2402      382851 :       av = avma; p1 = gen_0;
    2403     2713109 :       for (j=ex; j<=minss(i,mi); j++)
    2404     2330258 :         p1 = gadd(p1, gmulgs(gmul(gel(xd,j),gel(yd,i-j)),j));
    2405      382851 :       gel(yd,i) = gerepileupto(av, gdivgs(p1,i));
    2406             :     }
    2407       40005 :     return y;
    2408             :   }
    2409       16072 :   av = avma;
    2410       16072 :   return gerepileupto(av, gmul(gexp(gel(x,2),prec), serexp(serchop0(x),prec)));
    2411             : }
    2412             : 
    2413             : GEN
    2414     8582372 : gexp(GEN x, long prec)
    2415             : {
    2416     8582372 :   switch(typ(x))
    2417             :   {
    2418     6921353 :     case t_REAL: return mpexp(x);
    2419     1414313 :     case t_COMPLEX: return cxexp(x,prec);
    2420          56 :     case t_PADIC: return Qp_exp(x);
    2421             :     default:
    2422             :     {
    2423      246650 :       pari_sp av = avma;
    2424             :       GEN y;
    2425      246650 :       if (!(y = toser_i(x))) break;
    2426       58856 :       return gerepileupto(av, serexp(y,prec));
    2427             :     }
    2428             :   }
    2429      187794 :   return trans_eval("exp",gexp,x,prec);
    2430             : }
    2431             : 
    2432             : /********************************************************************/
    2433             : /**                                                                **/
    2434             : /**                           AGM(X, Y)                            **/
    2435             : /**                                                                **/
    2436             : /********************************************************************/
    2437             : static int
    2438     4751728 : agmr_gap(GEN a, GEN b, long L)
    2439             : {
    2440     4751728 :   GEN d = subrr(b, a);
    2441     4751728 :   return (signe(d) && expo(d) - expo(b) >= L);
    2442             : }
    2443             : /* assume x > 0 */
    2444             : static GEN
    2445      338072 : agm1r_abs(GEN x)
    2446             : {
    2447      338072 :   long l = realprec(x), L = 5-prec2nbits(l);
    2448      338072 :   GEN a1, b1, y = cgetr(l);
    2449      338072 :   pari_sp av = avma;
    2450             : 
    2451      338072 :   a1 = addrr(real_1(l), x); shiftr_inplace(a1, -1);
    2452      338072 :   b1 = sqrtr_abs(x);
    2453     5089800 :   while (agmr_gap(a1,b1,L))
    2454             :   {
    2455     4413656 :     GEN a = a1;
    2456     4413656 :     a1 = addrr(a,b1); shiftr_inplace(a1, -1);
    2457     4413656 :     b1 = sqrtr_abs(mulrr(a,b1));
    2458             :   }
    2459      338072 :   affrr_fixlg(a1,y); set_avma(av); return y;
    2460             : }
    2461             : 
    2462             : struct agmcx_gap_t { long L, ex, cnt; };
    2463             : 
    2464             : static void
    2465       14579 : agmcx_init(GEN x, long *prec, struct agmcx_gap_t *S)
    2466             : {
    2467       14579 :   long l = precision(x);
    2468       14579 :   if (l) *prec = l;
    2469       14579 :   S->L = 1-prec2nbits(*prec);
    2470       14579 :   S->cnt = 0;
    2471       14579 :   S->ex = LONG_MAX;
    2472       14579 : }
    2473             : 
    2474             : static long
    2475       14579 : agmcx_a_b(GEN x, GEN *a1, GEN *b1, long prec)
    2476             : {
    2477       14579 :   long rotate = 0;
    2478       14579 :   if (gsigne(real_i(x))<0)
    2479             :   { /* Rotate by +/-Pi/2, so that the choice of the principal square
    2480             :      * root gives the optimal AGM. So a1 = +/-I*a1, b1=sqrt(-x). */
    2481        1477 :     if (gsigne(imag_i(x))<0) { *a1=mulcxI(*a1);  rotate=-1; }
    2482         987 :     else                     { *a1=mulcxmI(*a1); rotate=1; }
    2483        1477 :     x = gneg(x);
    2484             :   }
    2485       14579 :   *b1 = gsqrt(x, prec);
    2486       14579 :   return rotate;
    2487             : }
    2488             : /* return 0 if we must stop the AGM loop (a=b or a ~ b), 1 otherwise */
    2489             : static int
    2490      207027 : agmcx_gap(GEN a, GEN b, struct agmcx_gap_t *S)
    2491             : {
    2492      207027 :   GEN d = gsub(b, a);
    2493      207027 :   long ex = S->ex;
    2494      207027 :   S->ex = gexpo(d);
    2495      207027 :   if (gequal0(d) || S->ex - gexpo(b) < S->L) return 0;
    2496             :   /* if (S->ex >= ex) we're no longer making progress; twice in a row */
    2497      195904 :   if (S->ex < ex) S->cnt = 0;
    2498             :   else
    2499        6961 :     if (S->cnt++) return 0;
    2500      192448 :   return 1;
    2501             : }
    2502             : static GEN
    2503       14537 : agm1cx(GEN x, long prec)
    2504             : {
    2505             :   struct agmcx_gap_t S;
    2506             :   GEN a1, b1;
    2507       14537 :   pari_sp av = avma;
    2508             :   long rotate;
    2509       14537 :   agmcx_init(x, &prec, &S);
    2510       14537 :   a1 = gtofp(gmul2n(gadd(real_1(prec), x), -1), prec);
    2511       14537 :   rotate = agmcx_a_b(x, &a1, &b1, prec);
    2512      221264 :   while (agmcx_gap(a1,b1,&S))
    2513             :   {
    2514      192190 :     GEN a = a1;
    2515      192190 :     a1 = gmul2n(gadd(a,b1),-1);
    2516      192190 :     b1 = gsqrt(gmul(a,b1), prec);
    2517             :   }
    2518       14537 :   if (rotate) a1 = rotate>0 ? mulcxI(a1):mulcxmI(a1);
    2519       14537 :   return gerepilecopy(av,a1);
    2520             : }
    2521             : 
    2522             : GEN
    2523          42 : zellagmcx(GEN a0, GEN b0, GEN r, GEN t, long prec)
    2524             : {
    2525             :   struct agmcx_gap_t S;
    2526          42 :   pari_sp av = avma;
    2527          42 :   GEN x = gdiv(a0, b0), a1, b1;
    2528             :   long rotate;
    2529          42 :   agmcx_init(x, &prec, &S);
    2530          42 :   a1 = gtofp(gmul2n(gadd(real_1(prec), x), -1), prec);
    2531          42 :   r = gsqrt(gdiv(gmul(a1,gaddgs(r, 1)),gadd(r, x)), prec);
    2532          42 :   t = gmul(r, t);
    2533          42 :   rotate = agmcx_a_b(x, &a1, &b1, prec);
    2534         342 :   while (agmcx_gap(a1,b1,&S))
    2535             :   {
    2536         258 :     GEN a = a1, b = b1;
    2537         258 :     a1 = gmul2n(gadd(a,b),-1);
    2538         258 :     b1 = gsqrt(gmul(a,b), prec);
    2539         258 :     r = gsqrt(gdiv(gmul(a1,gaddgs(r, 1)),gadd(gmul(b, r), a )), prec);
    2540         258 :     t = gmul(r, t);
    2541             :   }
    2542          42 :   if (rotate) a1 = rotate>0 ? mulcxI(a1):mulcxmI(a1);
    2543          42 :   a1 = gmul(a1, b0);
    2544          42 :   t = gatan(gdiv(a1,t), prec);
    2545             :   /* send t to the fundamental domain if necessary */
    2546          42 :   if (gsigne(real_i(t))<0) t = gadd(t, mppi(prec));
    2547          42 :   return gerepileupto(av,gdiv(t,a1));
    2548             : }
    2549             : 
    2550             : static long
    2551          49 : ser_cmp_expo(GEN A, GEN B)
    2552             : {
    2553          49 :   long e = -(long)HIGHEXPOBIT, d = valp(B) - valp(A);
    2554          49 :   long i, la = lg(A), v = varn(B);
    2555        9849 :   for (i = 2; i < la; i++)
    2556             :   {
    2557        9800 :     GEN a = gel(A,i), b;
    2558             :     long ei;
    2559        9800 :     if (isexactzero(a)) continue;
    2560        9800 :     b = polcoef_i(B, i-2 + d, v);
    2561        9800 :     ei = gexpo(a);
    2562        9800 :     if (!isexactzero(b)) ei -= gexpo(b);
    2563        9800 :     e = maxss(e, ei);
    2564             :   }
    2565          49 :   return e;
    2566             : }
    2567             : 
    2568             : static GEN
    2569          14 : ser_agm1(GEN y, long prec)
    2570             : {
    2571          14 :   GEN a1 = y, b1 = gen_1;
    2572          14 :   long l = lg(y)-2, l2 = 6-prec2nbits(prec), eold = LONG_MAX;
    2573             :   for(;;)
    2574          63 :   {
    2575          77 :     GEN a = a1, p1;
    2576          77 :     a1 = gmul2n(gadd(a,b1),-1);
    2577          77 :     b1 = gsqrt(gmul(a,b1), prec);
    2578          77 :     p1 = gsub(b1,a1);
    2579          77 :     if (isinexactreal(p1))
    2580             :     {
    2581          49 :       long e = ser_cmp_expo(p1, b1);
    2582          49 :       if (e < l2 || e >= eold) break;
    2583          42 :       eold = e;
    2584             :     }
    2585             :     else
    2586             :     {
    2587          28 :       long ep = valp(p1)-valp(b1);
    2588          28 :       if (ep >= l && gequal0(p1)) break;
    2589             :     }
    2590             :   }
    2591          14 :   return a1;
    2592             : }
    2593             : 
    2594             : /* agm(1,x) */
    2595             : static GEN
    2596       10768 : agm1(GEN x, long prec)
    2597             : {
    2598             :   GEN y;
    2599             :   pari_sp av;
    2600             : 
    2601       10768 :   if (gequal0(x)) return gcopy(x);
    2602       10768 :   switch(typ(x))
    2603             :   {
    2604             :     case t_INT:
    2605          28 :       if (!is_pm1(x)) break;
    2606          21 :       return (signe(x) > 0)? real_1(prec): real_0(prec);
    2607             : 
    2608        5245 :     case t_REAL: return signe(x) > 0? agm1r_abs(x): agm1cx(x, prec);
    2609             : 
    2610             :     case t_COMPLEX:
    2611        5362 :       if (gequal0(gel(x,2))) return agm1(gel(x,1), prec);
    2612        5173 :       return agm1cx(x, prec);
    2613             : 
    2614             :     case t_PADIC:
    2615             :     {
    2616          14 :       GEN a1 = x, b1 = gen_1;
    2617          14 :       long l = precp(x);
    2618          14 :       av = avma;
    2619             :       for(;;)
    2620          14 :       {
    2621          28 :         GEN a = a1, p1;
    2622             :         long ep;
    2623          28 :         a1 = gmul2n(gadd(a,b1),-1);
    2624          28 :         a = gmul(a,b1);
    2625          28 :         b1 = Qp_sqrt(a); if (!b1) pari_err_SQRTN("Qp_sqrt",a);
    2626          21 :         p1 = gsub(b1,a1); ep = valp(p1)-valp(b1);
    2627          21 :         if (ep<=0) { b1 = gneg_i(b1); p1 = gsub(b1,a1); ep=valp(p1)-valp(b1); }
    2628          21 :         if (ep >= l || gequal0(p1)) return gerepilecopy(av,a1);
    2629             :       }
    2630             :     }
    2631             : 
    2632             :     default:
    2633         119 :       av = avma; if (!(y = toser_i(x))) break;
    2634          14 :       return gerepilecopy(av, ser_agm1(y, prec));
    2635             :   }
    2636         112 :   return trans_eval("agm",agm1,x,prec);
    2637             : }
    2638             : 
    2639             : GEN
    2640       10467 : agm(GEN x, GEN y, long prec)
    2641             : {
    2642             :   pari_sp av;
    2643       10467 :   if (is_matvec_t(typ(y)))
    2644             :   {
    2645          14 :     if (is_matvec_t(typ(x))) pari_err_TYPE2("agm",x,y);
    2646           7 :     swap(x, y);
    2647             :   }
    2648       10460 :   if (gequal0(y)) return gcopy(y);
    2649       10460 :   av = avma;
    2650       10460 :   return gerepileupto(av, gmul(y, agm1(gdiv(x,y), prec)));
    2651             : }
    2652             : 
    2653             : /********************************************************************/
    2654             : /**                                                                **/
    2655             : /**                             LOG(X)                             **/
    2656             : /**                                                                **/
    2657             : /********************************************************************/
    2658             : /* atanh(u/v) using binary splitting */
    2659             : static GEN
    2660        5366 : atanhQ_split(ulong u, ulong v, long prec)
    2661             : {
    2662             :   long i, nmax;
    2663        5366 :   GEN u2 = sqru(u), v2 = sqru(v);
    2664        5366 :   double d = ((double)v) / u;
    2665             :   struct abpq_res R;
    2666             :   struct abpq A;
    2667             :   /* satisfies (2n+1) (v/u)^2n > 2^bitprec */
    2668        5366 :   nmax = (long)(prec2nbits(prec) / (2*log2(d)));
    2669        5366 :   abpq_init(&A, nmax);
    2670        5366 :   A.a[0] = A.b[0] = gen_1;
    2671        5366 :   A.p[0] = utoipos(u);
    2672        5366 :   A.q[0] = utoipos(v);
    2673      471547 :   for (i = 1; i <= nmax; i++)
    2674             :   {
    2675      466181 :     A.a[i] = gen_1;
    2676      466181 :     A.b[i] = utoipos((i<<1)+1);
    2677      466181 :     A.p[i] = u2;
    2678      466181 :     A.q[i] = v2;
    2679             :   }
    2680        5366 :   abpq_sum(&R, 0, nmax, &A);
    2681        5366 :   return rdivii(R.T, mulii(R.B,R.Q),prec);
    2682             : }
    2683             : /* log(2) = 10*atanh(1/17)+4*atanh(13/499); faster than logagmr_abs()
    2684             :  * and Pi/2M(1,4/2^n) ~ n log(2) */
    2685             : static GEN
    2686        2683 : log2_split(long prec)
    2687             : {
    2688        2683 :   GEN u = atanhQ_split(1, 17, prec);
    2689        2683 :   GEN v = atanhQ_split(13, 499, prec);
    2690        2683 :   shiftr_inplace(v, 2);
    2691        2683 :   return addrr(mulur(10, u), v);
    2692             : }
    2693             : GEN
    2694    13352038 : constlog2(long prec)
    2695             : {
    2696             :   pari_sp av;
    2697             :   GEN tmp;
    2698    13352038 :   if (glog2 && realprec(glog2) >= prec) return glog2;
    2699             : 
    2700        2683 :   tmp = cgetr_block(prec);
    2701        2683 :   av = avma;
    2702        2683 :   affrr(log2_split(prec+EXTRAPRECWORD), tmp);
    2703        2683 :   swap_clone(&glog2,tmp);
    2704        2683 :   set_avma(av); return glog2;
    2705             : }
    2706             : 
    2707             : GEN
    2708    13352038 : mplog2(long prec) { return rtor(constlog2(prec), prec); }
    2709             : 
    2710             : /* dont check that q != 2^expo(q), done in logr_abs */
    2711             : static GEN
    2712      332869 : logagmr_abs(GEN q)
    2713             : {
    2714      332869 :   long prec = realprec(q), e = expo(q), lim;
    2715      332869 :   GEN z = cgetr(prec), y, Q, _4ovQ;
    2716      332869 :   pari_sp av = avma;
    2717             : 
    2718      332869 :   incrprec(prec);
    2719      332869 :   lim = prec2nbits(prec) >> 1;
    2720      332869 :   Q = rtor(q,prec);
    2721      332869 :   shiftr_inplace(Q,lim-e); setsigne(Q,1);
    2722             : 
    2723      332869 :   _4ovQ = invr(Q); shiftr_inplace(_4ovQ, 2); /* 4/Q */
    2724             :   /* Pi / 2agm(1, 4/Q) ~ log(Q), q = Q * 2^(e-lim) */
    2725      332869 :   y = divrr(Pi2n(-1, prec), agm1r_abs(_4ovQ));
    2726      332869 :   y = addrr(y, mulsr(e - lim, mplog2(prec)));
    2727      332869 :   affrr_fixlg(y, z); set_avma(av); return z;
    2728             : }
    2729             : 
    2730             : /* sum_{k >= 0} y^(2k+1) / (2k+1), y close to 0 */
    2731             : static GEN
    2732     3211974 : logr_aux(GEN y)
    2733             : {
    2734     3211974 :   long k, L = realprec(y); /* should be ~ l+1 - (k-2) */
    2735             :   /* log(x) = log(1+y) - log(1-y) = 2 sum_{k odd} y^k / k
    2736             :    * Truncate the sum at k = 2n+1, the remainder is
    2737             :    *   2 sum_{k >= 2n+3} y^k / k < 2y^(2n+3) / (2n+3)(1-y) < y^(2n+3)
    2738             :    * We want y^(2n+3) < y 2^(-prec2nbits(L)), hence
    2739             :    *   n+1 > -prec2nbits(L) /-log_2(y^2) */
    2740     3211974 :   double d = -2*dbllog2r(y); /* ~ -log_2(y^2) */
    2741     3211974 :   k = (long)(2*(prec2nbits(L) / d));
    2742     3211974 :   k |= 1;
    2743     3211974 :   if (k >= 3)
    2744             :   {
    2745     3190466 :     GEN T, S = cgetr(L), y2 = sqrr(y), unr = real_1(L);
    2746     3190466 :     pari_sp av = avma;
    2747     3190466 :     long s = 0, incs = (long)d, l1 = nbits2prec((long)d);
    2748     3190466 :     setprec(S,  l1);
    2749     3190466 :     setprec(unr,l1); affrr(divru(unr,k), S);
    2750    51992282 :     for (k -= 2;; k -= 2) /* k = 2n+1, ..., 1 */
    2751             :     { /* S = y^(2n+1-k)/(2n+1) + ... + 1 / k */
    2752   100794098 :       setprec(y2, l1); T = mulrr(S,y2);
    2753    51992282 :       if (k == 1) break;
    2754             : 
    2755    48801816 :       l1 += dvmdsBIL(s + incs, &s); if (l1>L) l1=L;
    2756    48801816 :       setprec(S, l1);
    2757    48801816 :       setprec(unr,l1);
    2758    48801816 :       affrr(addrr(divru(unr, k), T), S); set_avma(av);
    2759             :     }
    2760             :     /* k = 1 special-cased for eficiency */
    2761     3190466 :     y = mulrr(y, addsr(1,T)); /* = log(X)/2 */
    2762             :   }
    2763     3211974 :   return y;
    2764             : }
    2765             : /*return log(|x|), assuming x != 0 */
    2766             : GEN
    2767     3772525 : logr_abs(GEN X)
    2768             : {
    2769     3772525 :   long EX, L, m, k, a, b, l = realprec(X);
    2770             :   GEN z, x, y;
    2771             :   ulong u;
    2772             :   double d;
    2773             : 
    2774             :  /* Assuming 1 < x < 2, we want delta = x-1, 1-x/2, 1-1/x, or 2/x-1 small.
    2775             :   * We have 2/x-1 > 1-x/2, 1-1/x < x-1. So one should be choosing between
    2776             :   * 1-1/x and 1-x/2 ( crossover sqrt(2), worse ~ 0.29 ). To avoid an inverse,
    2777             :   * we choose between x-1 and 1-x/2 ( crossover 4/3, worse ~ 0.33 ) */
    2778     3772525 :   EX = expo(X);
    2779     3772525 :   u = uel(X,2);
    2780     3772525 :   k = 2;
    2781     3772525 :   if (u > (~0UL / 3) * 2) { /* choose 1-x/2 */
    2782     2070795 :     EX++; u = ~u;
    2783     2070795 :     while (!u && ++k < l) { u = uel(X,k); u = ~u; }
    2784             :   } else { /* choose x - 1 */
    2785     1701730 :     u &= ~HIGHBIT; /* u - HIGHBIT, assuming HIGHBIT set */
    2786     1701730 :     while (!u && ++k < l) u = uel(X,k);
    2787             :   }
    2788     3772525 :   if (k == l) return EX? mulsr(EX, mplog2(l)): real_0(l);
    2789     3544801 :   a = prec2nbits(k) + bfffo(u); /* ~ -log2 |1-x| */
    2790     3544801 :   L = l+1;
    2791     3544801 :   b = prec2nbits(L - (k-2)); /* take loss of accuracy into account */
    2792     6050892 :   if (b > 24*a*log2(L) &&
    2793     2838960 :       prec2nbits(l) > prec2nbits(LOGAGM_LIMIT)) return logagmr_abs(X);
    2794             : 
    2795     3211932 :   z = cgetr(EX? l: l - (k-2));
    2796             : 
    2797             :  /* Multiplication is quadratic in this range (l is small, otherwise we
    2798             :   * use AGM). Set Y = x^(1/2^m), y = (Y - 1) / (Y + 1) and compute truncated
    2799             :   * series sum y^(2k+1)/(2k+1): the costs is less than
    2800             :   *    m b^2 + sum_{k <= n} ((2k+1) e + BITS_IN_LONG)^2
    2801             :   * bit operations with |x-1| <  2^(1-a), |Y| < 2^(1-e) , m = e-a and b bits of
    2802             :   * accuracy needed (+ BITS_IN_LONG since bit accuracies increase by
    2803             :   * increments of BITS_IN_LONG), so
    2804             :   * 4n^3/3 e^2 + n^2 2e BITS_IN_LONG+ n BITS_IN_LONG ~ m b^2, with n ~ b/2e
    2805             :   * or b/6e + BITS_IN_LONG/2e + BITS_IN_LONG/2be ~ m
    2806             :   *    B := (b / 6 + BITS_IN_LONG/2 + BITS_IN_LONG^2 / 2b) ~ m(m+a)
    2807             :   *     m = min( -a/2 + sqrt(a^2/4 + B),  b - a )
    2808             :   * NB: e ~ (b/6)^(1/2) as b -> oo
    2809             :   * Instead of the above pessimistic estimate for the cost of the sum, use
    2810             :   * optimistic estimate (BITS_IN_LONG -> 0) */
    2811     3211932 :   d = -a/2.; m = (long)(d + sqrt(d*d + b/6)); /* >= 0 */
    2812             : 
    2813     3211932 :   if (m > b-a) m = b-a;
    2814     3211932 :   if (m < 0.2*a) m = 0; else L += nbits2extraprec(m);
    2815     3211932 :   x = rtor(X,L);
    2816     3211932 :   setsigne(x,1); shiftr_inplace(x,-EX);
    2817             :   /* 2/3 < x < 4/3 */
    2818     3211932 :   for (k=1; k<=m; k++) x = sqrtr_abs(x);
    2819             : 
    2820     3211932 :   y = divrr(subrs(x,1), addrs(x,1)); /* = (x-1) / (x+1), close to 0 */
    2821     3211932 :   y = logr_aux(y); /* log(1+y) - log(1-y) = log(x) */
    2822     3211932 :   shiftr_inplace(y, m + 1);
    2823     3211932 :   if (EX) y = addrr(y, mulsr(EX, mplog2(l+1)));
    2824     3211932 :   affrr_fixlg(y, z); set_avma((pari_sp)z); return z;
    2825             : }
    2826             : 
    2827             : /* assume Im(q) != 0 and precision(q) >= prec. Compute log(q) with accuracy
    2828             :  * prec [disregard input accuracy] */
    2829             : GEN
    2830        9322 : logagmcx(GEN q, long prec)
    2831             : {
    2832        9322 :   GEN z = cgetc(prec), y, Q, a, b;
    2833             :   long lim, e, ea, eb;
    2834        9322 :   pari_sp av = avma;
    2835        9322 :   int neg = 0;
    2836             : 
    2837        9322 :   incrprec(prec);
    2838        9322 :   if (gsigne(gel(q,1)) < 0) { q = gneg(q); neg = 1; }
    2839        9322 :   lim = prec2nbits(prec) >> 1;
    2840        9322 :   Q = gtofp(q, prec);
    2841        9322 :   a = gel(Q,1);
    2842        9322 :   b = gel(Q,2);
    2843        9322 :   if (gequal0(a)) {
    2844           0 :     affrr_fixlg(logr_abs(b), gel(z,1));
    2845           0 :     y = Pi2n(-1, prec);
    2846           0 :     if (signe(b) < 0) setsigne(y, -1);
    2847           0 :     affrr_fixlg(y, gel(z,2)); set_avma(av); return z;
    2848             :   }
    2849        9322 :   ea = expo(a);
    2850        9322 :   eb = expo(b);
    2851        9322 :   e = ea <= eb ? lim - eb : lim - ea;
    2852        9322 :   shiftr_inplace(a, e);
    2853        9322 :   shiftr_inplace(b, e);
    2854             : 
    2855             :   /* Pi / 2agm(1, 4/Q) ~ log(Q), q = Q * 2^e */
    2856        9322 :   y = gdiv(Pi2n(-1, prec), agm1cx( gdivsg(4, Q), prec ));
    2857        9322 :   a = gel(y,1);
    2858        9322 :   b = gel(y,2);
    2859        9322 :   a = addrr(a, mulsr(-e, mplog2(prec)));
    2860        9322 :   if (realprec(a) <= LOWDEFAULTPREC) a = real_0_bit(expo(a));
    2861       15754 :   if (neg) b = gsigne(b) <= 0? gadd(b, mppi(prec))
    2862        6432 :                              : gsub(b, mppi(prec));
    2863        9322 :   affrr_fixlg(a, gel(z,1));
    2864        9322 :   affrr_fixlg(b, gel(z,2)); set_avma(av); return z;
    2865             : }
    2866             : 
    2867             : GEN
    2868      107229 : mplog(GEN x)
    2869             : {
    2870      107229 :   if (signe(x)<=0) pari_err_DOMAIN("mplog", "argument", "<=", gen_0, x);
    2871      107229 :   return logr_abs(x);
    2872             : }
    2873             : 
    2874             : /* pe = p^e, p prime, 0 < x < pe a t_INT coprime to p. Return the (p-1)-th
    2875             :  * root of 1 in (Z/pe)^* congruent to x mod p, resp x mod 4 if p = 2.
    2876             :  * Simplified form of Zp_sqrtnlift: 1/(p-1) is trivial to compute */
    2877             : GEN
    2878        1344 : Zp_teichmuller(GEN x, GEN p, long e, GEN pe)
    2879             : {
    2880             :   GEN q, z, p1;
    2881             :   pari_sp av;
    2882             :   ulong mask;
    2883        1344 :   if (absequaliu(p,2)) return (mod4(x) & 2)? subiu(pe,1): gen_1;
    2884        1344 :   if (e == 1) return icopy(x);
    2885        1344 :   av = avma;
    2886        1344 :   p1 = subiu(p, 1);
    2887        1344 :   mask = quadratic_prec_mask(e);
    2888        1344 :   q = p; z = remii(x, p);
    2889        7868 :   while (mask > 1)
    2890             :   { /* Newton iteration solving z^{1 - p} = 1, z = x (mod p) */
    2891        5180 :     GEN w, t, qold = q;
    2892        5180 :     if (mask <= 3) /* last iteration */
    2893        1344 :       q = pe;
    2894             :     else
    2895             :     {
    2896        3836 :       q = sqri(q);
    2897        3836 :       if (mask & 1) q = diviiexact(q, p);
    2898             :     }
    2899        5180 :     mask >>= 1;
    2900             :     /* q <= qold^2 */
    2901        5180 :     if (lgefint(q) == 3)
    2902             :     {
    2903        5032 :       ulong Z = uel(z,2), Q = uel(q,2), P1 = uel(p1,2);
    2904        5032 :       ulong W = (Q-1) / P1; /* -1/(p-1) + O(qold) */
    2905        5032 :       ulong T = Fl_mul(W, Fl_powu(Z,P1,Q) - 1, Q);
    2906        5032 :       Z = Fl_mul(Z, 1 + T, Q);
    2907        5032 :       z = utoi(Z);
    2908             :     }
    2909             :     else
    2910             :     {
    2911         148 :       w = diviiexact(subiu(qold,1),p1); /* -1/(p-1) + O(qold) */
    2912         148 :       t = Fp_mul(w, subiu(Fp_pow(z,p1,q), 1), q);
    2913         148 :       z = Fp_mul(z, addui(1,t), q);
    2914             :     }
    2915             :   }
    2916        1344 :   return gerepileuptoint(av, z);
    2917             : }
    2918             : 
    2919             : GEN
    2920        1225 : teichmullerinit(long p, long n)
    2921             : {
    2922             :   GEN t, pn, g, v;
    2923             :   ulong gp, tp;
    2924             :   long a, m;
    2925             : 
    2926        1225 :   if (p == 2) return mkvec(gen_1);
    2927        1225 :   if (!uisprime(p)) pari_err_PRIME("teichmullerinit",utoipos(p));
    2928             : 
    2929        1225 :   m = p >> 1; /* (p-1)/2 */
    2930        1225 :   tp= gp= pgener_Fl(p); /* order (p-1), gp^m = -1 */
    2931        1225 :   pn = powuu(p, n);
    2932        1225 :   v = cgetg(p, t_VEC);
    2933        1225 :   t = g = Zp_teichmuller(utoipos(gp), utoipos(p), n, pn);
    2934        1225 :   gel(v, 1) = gen_1;
    2935        1225 :   gel(v, p-1) = subiu(pn,1);
    2936        3031 :   for (a = 1; a < m; a++)
    2937             :   {
    2938        1806 :     gel(v, tp) = t;
    2939        1806 :     gel(v, p - tp) = Fp_neg(t, pn); /* g^(m+a) = -g^a */
    2940        1806 :     if (a < m-1)
    2941             :     {
    2942        1029 :       t = Fp_mul(t, g, pn); /* g^(a+1) */
    2943        1029 :       tp = Fl_mul(tp, gp, p); /* t mod p  */
    2944             :     }
    2945             :   }
    2946        1225 :   return v;
    2947             : }
    2948             : 
    2949             : /* tab from teichmullerinit or NULL */
    2950             : GEN
    2951         238 : teichmuller(GEN x, GEN tab)
    2952             : {
    2953             :   GEN p, q, y, z;
    2954         238 :   long n, tx = typ(x);
    2955             : 
    2956         238 :   if (!tab)
    2957             :   {
    2958         126 :     if (tx == t_VEC && lg(x) == 3)
    2959             :     {
    2960           7 :       p = gel(x,1);
    2961           7 :       q = gel(x,2);
    2962           7 :       if (typ(p) == t_INT && typ(q) == t_INT)
    2963           7 :         return teichmullerinit(itos(p), itos(q));
    2964             :     }
    2965             :   }
    2966         112 :   else if (typ(tab) != t_VEC) pari_err_TYPE("teichmuller",tab);
    2967         231 :   if (tx!=t_PADIC) pari_err_TYPE("teichmuller",x);
    2968         231 :   z = gel(x,4);
    2969         231 :   if (!signe(z)) return gcopy(x);
    2970         231 :   p = gel(x,2);
    2971         231 :   q = gel(x,3);
    2972         231 :   n = precp(x);
    2973         231 :   y = cgetg(5,t_PADIC);
    2974         231 :   y[1] = evalprecp(n) | _evalvalp(0);
    2975         231 :   gel(y,2) = icopy(p);
    2976         231 :   gel(y,3) = icopy(q);
    2977         231 :   if (tab)
    2978             :   {
    2979         112 :     ulong pp = itou_or_0(p);
    2980         112 :     if (lg(tab) != (long)pp) pari_err_TYPE("teichmuller",tab);
    2981         112 :     z = gel(tab, umodiu(z, pp));
    2982         112 :     if (typ(z) != t_INT) pari_err_TYPE("teichmuller",tab);
    2983         112 :     z = remii(z, q);
    2984             :   }
    2985             :   else
    2986         119 :     z = Zp_teichmuller(z, p, n, q);
    2987         231 :   gel(y,4) = z;
    2988         231 :   return y;
    2989             : }
    2990             : GEN
    2991           0 : teich(GEN x) { return teichmuller(x, NULL); }
    2992             : 
    2993             : GEN
    2994     3052398 : glog(GEN x, long prec)
    2995             : {
    2996             :   pari_sp av, tetpil;
    2997             :   GEN y, p1;
    2998             :   long l;
    2999             : 
    3000     3052398 :   switch(typ(x))
    3001             :   {
    3002             :     case t_REAL:
    3003     2036667 :       if (signe(x) >= 0)
    3004             :       {
    3005     1761367 :         if (!signe(x)) pari_err_DOMAIN("log", "argument", "=", gen_0, x);
    3006     1761353 :         return logr_abs(x);
    3007             :       }
    3008      275300 :       retmkcomplex(logr_abs(x), mppi(realprec(x)));
    3009             : 
    3010             :     case t_FRAC:
    3011             :     {
    3012             :       GEN a, b;
    3013             :       long e1, e2;
    3014      127675 :       av = avma;
    3015      127675 :       a = gel(x,1);
    3016      127675 :       b = gel(x,2);
    3017      127675 :       e1 = expi(subii(a,b)); e2 = expi(b);
    3018      127675 :       if (e2 > e1) prec += nbits2nlong(e2 - e1);
    3019      127675 :       x = fractor(x, prec);
    3020      127675 :       return gerepileupto(av, glog(x, prec));
    3021             :     }
    3022             :     case t_COMPLEX:
    3023      592099 :       if (ismpzero(gel(x,2))) return glog(gel(x,1), prec);
    3024      587538 :       l = precision(x); if (l > prec) prec = l;
    3025      587538 :       if (ismpzero(gel(x,1)))
    3026             :       {
    3027        5367 :         GEN a = gel(x,2), b;
    3028        5367 :         av = avma; b = Pi2n(-1,prec);
    3029        5367 :         if (gsigne(a) < 0) { setsigne(b, -1); a = gabs(a,prec); }
    3030        5367 :         a = isint1(a) ? gen_0: glog(a,prec);
    3031        5367 :         return gerepilecopy(av, mkcomplex(a, b));
    3032             :       }
    3033      582171 :       if (prec >= LOGAGMCX_LIMIT) return logagmcx(x, prec);
    3034      573045 :       y = cgetg(3,t_COMPLEX);
    3035      573045 :       gel(y,2) = garg(x,prec);
    3036      573045 :       av = avma; p1 = glog(cxnorm(x),prec); tetpil = avma;
    3037      573045 :       gel(y,1) = gerepile(av,tetpil,gmul2n(p1,-1)); return y;
    3038             : 
    3039         287 :     case t_PADIC: return Qp_log(x);
    3040             :     default:
    3041      295670 :       av = avma; if (!(y = toser_i(x))) break;
    3042         182 :       if (!signe(y)) pari_err_DOMAIN("log", "argument", "=", gen_0, x);
    3043         182 :       if (valp(y)) pari_err_DOMAIN("log", "series valuation", "!=", gen_0, x);
    3044         175 :       p1 = integser(gdiv(derivser(y), y)); /* log(y)' = y'/y */
    3045         175 :       if (!gequal1(gel(y,2))) p1 = gadd(p1, glog(gel(y,2),prec));
    3046         175 :       return gerepileupto(av, p1);
    3047             :   }
    3048      295488 :   return trans_eval("log",glog,x,prec);
    3049             : }
    3050             : 
    3051             : static GEN
    3052          42 : mplog1p(GEN x)
    3053             : {
    3054             :   long ex, a, b, l, L;
    3055          42 :   if (!signe(x)) return rcopy(x);
    3056          42 :   ex = expo(x); if (ex >= 0) return glog(addrs(x,1), 0);
    3057          42 :   a = -ex;
    3058          42 :   b = realprec(x); L = b+1;
    3059          42 :   if (b > a*log2(L) && prec2nbits(b) > prec2nbits(LOGAGM_LIMIT))
    3060             :   {
    3061           0 :     x = addrs(x,1); l = b + nbits2extraprec(a);
    3062           0 :     if (realprec(x) < l) x = rtor(x,l);
    3063           0 :     return logagmr_abs(x);
    3064             :   }
    3065          42 :   x = rtor(x, L);
    3066          42 :   x = logr_aux(divrr(x, addrs(x,2)));
    3067          42 :   if (realprec(x) > b) fixlg(x, b);
    3068          42 :   shiftr_inplace(x,1); return x;
    3069             : }
    3070             : 
    3071             : static GEN log1p_i(GEN x, long prec);
    3072             : static GEN
    3073          14 : cxlog1p(GEN x, long prec)
    3074             : {
    3075             :   pari_sp av;
    3076          14 :   GEN z, a, b = gel(x,2);
    3077             :   long l;
    3078          14 :   if (ismpzero(b)) return log1p_i(gel(x,1), prec);
    3079          14 :   l = precision(x); if (l > prec) prec = l;
    3080          14 :   if (prec >= LOGAGMCX_LIMIT) return logagmcx(gaddgs(x,1), prec);
    3081          14 :   a = gel(x,1);
    3082          14 :   z = cgetg(3,t_COMPLEX); av = avma;
    3083          14 :   a = gadd(gadd(gmul2n(a,1), gsqr(a)), gsqr(b));
    3084          14 :   a = log1p_i(a, prec); shiftr_inplace(a,-1);
    3085          14 :   gel(z,1) = gerepileupto(av, a);
    3086          14 :   gel(z,2) = garg(gaddgs(x,1),prec); return z;
    3087             : }
    3088             : static GEN
    3089          91 : log1p_i(GEN x, long prec)
    3090             : {
    3091          91 :   switch(typ(x))
    3092             :   {
    3093          42 :     case t_REAL: return mplog1p(x);
    3094          14 :     case t_COMPLEX: return cxlog1p(x, prec);
    3095           7 :     case t_PADIC: return Qp_log(gaddgs(x,1));
    3096             :     default:
    3097             :     {
    3098             :       long ey;
    3099             :       GEN y;
    3100          28 :       if (!(y = toser_i(x))) break;
    3101          21 :       ey = valp(y);
    3102          21 :       if (ey < 0) pari_err_DOMAIN("log1p","valuation", "<", gen_0, x);
    3103          21 :       if (gequal0(y)) return gcopy(y);
    3104          14 :       if (ey)
    3105           7 :         return glog(gaddgs(y,1),prec);
    3106             :       else
    3107             :       {
    3108           7 :         GEN a = gel(y,2), a1 = gaddgs(a,1);
    3109           7 :         y = gdiv(y, a1); gel(y,2) = gen_1;
    3110           7 :         return gadd(glog1p(a,prec), glog(y, prec));
    3111             :       }
    3112             :     }
    3113             :   }
    3114           7 :   return trans_eval("log1p",glog1p,x,prec);
    3115             : }
    3116             : GEN
    3117          77 : glog1p(GEN x, long prec)
    3118             : {
    3119          77 :   pari_sp av = avma;
    3120          77 :   return gerepileupto(av, log1p_i(x, prec));
    3121             : }
    3122             : /********************************************************************/
    3123             : /**                                                                **/
    3124             : /**                        SINE, COSINE                            **/
    3125             : /**                                                                **/
    3126             : /********************************************************************/
    3127             : 
    3128             : /* Reduce x0 mod Pi/2 to x in [-Pi/4, Pi/4]. Return cos(x)-1 */
    3129             : static GEN
    3130     4538951 : mpcosm1(GEN x, long *ptmod8)
    3131             : {
    3132     4538951 :   long a = expo(x), l = realprec(x), b, L, i, n, m, B;
    3133             :   GEN y, p2, x2;
    3134             :   double d;
    3135             : 
    3136     4538951 :   n = 0;
    3137     4538951 :   if (a >= 0)
    3138             :   {
    3139             :     long p;
    3140             :     GEN q;
    3141     3575577 :     if (a > 30)
    3142             :     {
    3143          92 :       GEN z, pitemp = Pi2n(-2, nbits2prec(a + 32));
    3144          92 :       z = addrr(x,pitemp); /* = x + Pi/4 */
    3145          92 :       if (expo(z) >= bit_prec(z) + 3) pari_err_PREC("mpcosm1");
    3146          92 :       shiftr_inplace(pitemp, 1);
    3147          92 :       q = floorr( divrr(z,pitemp) ); /* round ( x / (Pi/2) ) */
    3148          92 :       p = l+EXTRAPRECWORD; x = rtor(x,p);
    3149             :     } else {
    3150     3575485 :       q = stoi((long)floor(rtodbl(x) / (M_PI/2) + 0.5));
    3151     3575485 :       p = l;
    3152             :     }
    3153     3575577 :     if (signe(q))
    3154             :     {
    3155     3575577 :       x = subrr(x, mulir(q, Pi2n(-1,p))); /* x mod Pi/2  */
    3156     3575577 :       a = expo(x);
    3157     3575577 :       if (!signe(x) && a >= 0) pari_err_PREC("mpcosm1");
    3158     3575577 :       n = mod4(q); if (n && signe(q) < 0) n = 4 - n;
    3159             :     }
    3160             :   }
    3161             :   /* a < 0 */
    3162     4538951 :   b = signe(x); *ptmod8 = (b < 0)? 4 + n: n;
    3163     4538951 :   if (!b) return real_0_bit(expo(x)*2 - 1);
    3164             : 
    3165     4355673 :   b = prec2nbits(l);
    3166     4355673 :   if (b + 2*a <= 0) {
    3167      198674 :     y = sqrr(x); shiftr_inplace(y, -1); setsigne(y, -1);
    3168      198674 :     return y;
    3169             :   }
    3170             : 
    3171     4156999 :   y = cgetr(l);
    3172     4156999 :   B = b/6 + BITS_IN_LONG + (BITS_IN_LONG*BITS_IN_LONG/2)/ b;
    3173     4156999 :   d = a/2.; m = (long)(d + sqrt(d*d + B)); /* >= 0 ,*/
    3174     4156999 :   if (m < (-a) * 0.1) m = 0; /* not worth it */
    3175     4156999 :   L = l + nbits2extraprec(m);
    3176             : 
    3177     4156999 :   b += m;
    3178     4156999 :   d = 2.0 * (m-dbllog2r(x)-1/M_LN2); /* ~ 2( - log_2 Y - 1/log(2) ) */
    3179     4156999 :   n = (long)(b / d);
    3180     4156999 :   if (n > 1)
    3181     4125277 :     n = (long)(b / (d + log2((double)n+1))); /* log~constant in small ranges */
    3182     4156999 :   while (n*(d+log2((double)n+1)) < b) n++; /* expect few corrections */
    3183             : 
    3184             :  /* Multiplication is quadratic in this range (l is small, otherwise we
    3185             :   * use logAGM + Newton). Set Y = 2^(-e-a) x, compute truncated series
    3186             :   * sum Y^2k/(2k)!: this costs roughly
    3187             :   *   m b^2 + sum_{k <= n} (2k e + BITS_IN_LONG)^2
    3188             :   *   ~ floor(b/2e) b^2 / 3  + m b^2
    3189             :   * bit operations with |x| <  2^(1+a), |Y| < 2^(1-e) , m = e+a and b bits of
    3190             :   * accuracy needed, so
    3191             :   *    B := ( b / 6 + BITS_IN_LONG + BITS_IN_LONG^2 / 2b) ~ m(m-a)
    3192             :   * we want b ~ 6 m (m-a) or m~b+a hence
    3193             :   *     m = min( a/2 + sqrt(a^2/4 + b/6),  b/2 + a )
    3194             :   * NB1: e ~ (b/6)^(1/2) or b/2.
    3195             :   * NB2: We use b/4 instead of b/6 in the formula above: hand-optimized...
    3196             :   *
    3197             :   * Truncate the sum at k = n (>= 1), the remainder is
    3198             :   * < sum_{k >= n+1} Y^2k / 2k! < Y^(2n+2) / (2n+2)!(1-Y^2) < Y^(2n+2)/(2n+1)!
    3199             :   * We want ... <= Y^2 2^-b, hence -2n log_2 |Y| + log_2 (2n+1)! >= b
    3200             :   *   log n! ~ (n + 1/2) log(n+1) - (n+1) + log(2Pi)/2,
    3201             :   * error bounded by 1/6(n+1) <= 1/12. Finally, we want
    3202             :   * 2n (-1/log(2) - log_2 |Y| + log_2(2n+2)) >= b  */
    3203     4156999 :   x = rtor(x, L); shiftr_inplace(x, -m); setsigne(x, 1);
    3204     4156999 :   x2 = sqrr(x);
    3205     4156999 :   if (n == 1) { p2 = x2; shiftr_inplace(p2, -1); setsigne(p2, -1); } /*-Y^2/2*/
    3206             :   else
    3207             :   {
    3208     4156999 :     GEN unr = real_1(L);
    3209             :     pari_sp av;
    3210     4156999 :     long s = 0, l1 = nbits2prec((long)(d + n + 16));
    3211             : 
    3212     4156999 :     p2 = cgetr(L); av = avma;
    3213    26648589 :     for (i=n; i>=2; i--)
    3214             :     {
    3215             :       GEN p1;
    3216    22491590 :       setprec(x2,l1); p1 = divrunu(x2, 2*i-1);
    3217    22491590 :       l1 += dvmdsBIL(s - expo(p1), &s); if (l1>L) l1=L;
    3218    22491590 :       if (i != n) p1 = mulrr(p1,p2);
    3219    22491590 :       setprec(unr,l1); p1 = addrr_sign(unr,1, p1,-signe(p1));
    3220    22491590 :       setprec(p2,l1); affrr(p1,p2); set_avma(av);
    3221             :     }
    3222     4156999 :     shiftr_inplace(p2, -1); togglesign(p2); /* p2 := -p2/2 */
    3223     4156999 :     setprec(x2,L); p2 = mulrr(x2,p2);
    3224             :   }
    3225             :   /* Now p2 = sum {1<= i <=n} (-1)^i x^(2i) / (2i)! ~ cos(x) - 1 */
    3226    37925257 :   for (i=1; i<=m; i++)
    3227             :   { /* p2 = cos(x)-1 --> cos(2x)-1 */
    3228    33768258 :     p2 = mulrr(p2, addsr(2,p2));
    3229    33768258 :     shiftr_inplace(p2, 1);
    3230    33768258 :     if ((i & 31) == 0) p2 = gerepileuptoleaf((pari_sp)y, p2);
    3231             :   }
    3232     4156999 :   affrr_fixlg(p2,y); return y;
    3233             : }
    3234             : 
    3235             : /* sqrt (|1 - (1+x)^2|) = sqrt(|x*(x+2)|). Sends cos(x)-1 to |sin(x)| */
    3236             : static GEN
    3237     3067042 : mpaut(GEN x)
    3238             : {
    3239     3067042 :   pari_sp av = avma;
    3240     3067042 :   GEN t = mulrr(x, addsr(2,x)); /* != 0 */
    3241     3067042 :   if (!signe(t)) return real_0_bit(expo(t) >> 1);
    3242     2888426 :   return gerepileuptoleaf(av, sqrtr_abs(t));
    3243             : }
    3244             : 
    3245             : /********************************************************************/
    3246             : /**                            COSINE                              **/
    3247             : /********************************************************************/
    3248             : 
    3249             : GEN
    3250     2702158 : mpcos(GEN x)
    3251             : {
    3252             :   long mod8;
    3253             :   pari_sp av;
    3254             :   GEN y,p1;
    3255             : 
    3256     2702158 :   if (!signe(x)) {
    3257           8 :     long l = nbits2prec(-expo(x));
    3258           8 :     if (l < LOWDEFAULTPREC) l = LOWDEFAULTPREC;
    3259           8 :     return real_1(l);
    3260             :   }
    3261             : 
    3262     2702150 :   av = avma; p1 = mpcosm1(x,&mod8);
    3263     2702150 :   switch(mod8)
    3264             :   {
    3265      748978 :     case 0: case 4: y = addsr(1,p1); break;
    3266      677991 :     case 1: case 7: y = mpaut(p1); togglesign(y); break;
    3267      671600 :     case 2: case 6: y = subsr(-1,p1); break;
    3268      603581 :     default:        y = mpaut(p1); break; /* case 3: case 5: */
    3269             :   }
    3270     2702150 :   return gerepileuptoleaf(av, y);
    3271             : }
    3272             : 
    3273             : /* convert INT or FRAC to REAL, which is later reduced mod 2Pi : avoid
    3274             :  * cancellation */
    3275             : static GEN
    3276        6468 : tofp_safe(GEN x, long prec)
    3277             : {
    3278       11256 :   return (typ(x) == t_INT || gexpo(x) > 0)? gadd(x, real_0(prec))
    3279        8932 :                                           : fractor(x, prec);
    3280             : }
    3281             : 
    3282             : GEN
    3283      111959 : gcos(GEN x, long prec)
    3284             : {
    3285             :   pari_sp av;
    3286             :   GEN a, b, u, v, y, u1, v1;
    3287             :   long i;
    3288             : 
    3289      111959 :   switch(typ(x))
    3290             :   {
    3291      111070 :     case t_REAL: return mpcos(x);
    3292             :     case t_COMPLEX:
    3293          49 :       a = gel(x,1);
    3294          49 :       b = gel(x,2);
    3295          49 :       if (isintzero(a)) return gcosh(b, prec);
    3296          35 :       i = precision(x); if (i) prec = i;
    3297          35 :       y = cgetc(prec); av = avma;
    3298          35 :       if (typ(b) != t_REAL) b = gtofp(b, prec);
    3299          35 :       mpsinhcosh(b, &u1, &v1); u1 = mpneg(u1);
    3300          35 :       if (typ(a) != t_REAL) a = gtofp(a, prec);
    3301          35 :       mpsincos(a, &u, &v);
    3302          35 :       affrr_fixlg(gmul(v1,v), gel(y,1));
    3303          35 :       affrr_fixlg(gmul(u1,u), gel(y,2)); set_avma(av); return y;
    3304             : 
    3305             :     case t_INT: case t_FRAC:
    3306         756 :       y = cgetr(prec); av = avma;
    3307         756 :       affrr_fixlg(mpcos(tofp_safe(x,prec)), y); set_avma(av); return y;
    3308             : 
    3309          49 :     case t_PADIC: y = cos_p(x);
    3310          49 :       if (!y) pari_err_DOMAIN("gcos(t_PADIC)","argument","",gen_0,x);
    3311          42 :       return y;
    3312             : 
    3313             :     default:
    3314          35 :       av = avma; if (!(y = toser_i(x))) break;
    3315          28 :       if (gequal0(y)) return gerepileupto(av, gaddsg(1,y));
    3316          28 :       if (valp(y) < 0)
    3317           7 :         pari_err_DOMAIN("cos","valuation", "<", gen_0, x);
    3318          21 :       gsincos(y,&u,&v,prec);
    3319          21 :       return gerepilecopy(av,v);
    3320             :   }
    3321           7 :   return trans_eval("cos",gcos,x,prec);
    3322             : }
    3323             : /********************************************************************/
    3324             : /**                             SINE                               **/
    3325             : /********************************************************************/
    3326             : 
    3327             : GEN
    3328      157439 : mpsin(GEN x)
    3329             : {
    3330             :   long mod8;
    3331             :   pari_sp av;
    3332             :   GEN y,p1;
    3333             : 
    3334      157439 :   if (!signe(x)) return real_0_bit(expo(x));
    3335             : 
    3336      157283 :   av = avma; p1 = mpcosm1(x,&mod8);
    3337      157283 :   switch(mod8)
    3338             :   {
    3339       75978 :     case 0: case 6: y=mpaut(p1); break;
    3340       27779 :     case 1: case 5: y=addsr(1,p1); break;
    3341       29974 :     case 2: case 4: y=mpaut(p1); togglesign(y); break;
    3342       23552 :     default:        y=subsr(-1,p1); break; /* case 3: case 7: */
    3343             :   }
    3344      157283 :   return gerepileuptoleaf(av, y);
    3345             : }
    3346             : 
    3347             : GEN
    3348      188512 : gsin(GEN x, long prec)
    3349             : {
    3350             :   pari_sp av;
    3351             :   GEN a, b, u, v, y, v1, u1;
    3352             :   long i;
    3353             : 
    3354      188512 :   switch(typ(x))
    3355             :   {
    3356      152490 :     case t_REAL: return mpsin(x);
    3357             :     case t_COMPLEX:
    3358       31045 :       a = gel(x,1);
    3359       31045 :       b = gel(x,2);
    3360       31045 :       if (isintzero(a)) retmkcomplex(gen_0,gsinh(b,prec));
    3361       15722 :       i = precision(x); if (i) prec = i;
    3362       15722 :       y = cgetc(prec); av = avma;
    3363       15722 :       if (typ(b) != t_REAL) b = gtofp(b, prec);
    3364       15722 :       mpsinhcosh(b, &u1, &v1);
    3365       15722 :       if (typ(a) != t_REAL) a = gtofp(a, prec);
    3366       15722 :       mpsincos(a, &u, &v);
    3367       15722 :       affrr_fixlg(gmul(v1,u), gel(y,1));
    3368       15722 :       affrr_fixlg(gmul(u1,v), gel(y,2)); set_avma(av); return y;
    3369             : 
    3370             :     case t_INT: case t_FRAC:
    3371        4718 :       y = cgetr(prec); av = avma;
    3372        4718 :       affrr_fixlg(mpsin(tofp_safe(x,prec)), y); set_avma(av); return y;
    3373             : 
    3374          49 :     case t_PADIC: y = sin_p(x);
    3375          49 :       if (!y) pari_err_DOMAIN("gsin(t_PADIC)","argument","",gen_0,x);
    3376          42 :       return y;
    3377             : 
    3378             :     default:
    3379         210 :       av = avma; if (!(y = toser_i(x))) break;
    3380         203 :       if (gequal0(y)) return gerepilecopy(av, y);
    3381         203 :       if (valp(y) < 0)
    3382           7 :         pari_err_DOMAIN("sin","valuation", "<", gen_0, x);
    3383         196 :       gsincos(y,&u,&v,prec);
    3384         196 :       return gerepilecopy(av,u);
    3385             :   }
    3386           7 :   return trans_eval("sin",gsin,x,prec);
    3387             : }
    3388             : /********************************************************************/
    3389             : /**                       SINE, COSINE together                    **/
    3390             : /********************************************************************/
    3391             : 
    3392             : void
    3393     1678092 : mpsincos(GEN x, GEN *s, GEN *c)
    3394             : {
    3395             :   long mod8;
    3396             :   pari_sp av, tetpil;
    3397             :   GEN p1, *gptr[2];
    3398             : 
    3399     1678092 :   if (!signe(x))
    3400             :   {
    3401        3484 :     long e = expo(x);
    3402        3484 :     *s = real_0_bit(e);
    3403        3484 :     *c = e >= 0? real_0_bit(e): real_1_bit(-e);
    3404        6968 :     return;
    3405             :   }
    3406             : 
    3407     1674608 :   av=avma; p1=mpcosm1(x,&mod8); tetpil=avma;
    3408     1674608 :   switch(mod8)
    3409             :   {
    3410      412951 :     case 0: *c=addsr( 1,p1); *s=mpaut(p1); break;
    3411      123259 :     case 1: *s=addsr( 1,p1); *c=mpaut(p1); togglesign(*c); break;
    3412      353479 :     case 2: *c=subsr(-1,p1); *s=mpaut(p1); togglesign(*s); break;
    3413      113108 :     case 3: *s=subsr(-1,p1); *c=mpaut(p1); break;
    3414      234855 :     case 4: *c=addsr( 1,p1); *s=mpaut(p1); togglesign(*s); break;
    3415      140190 :     case 5: *s=addsr( 1,p1); *c=mpaut(p1); break;
    3416      179819 :     case 6: *c=subsr(-1,p1); *s=mpaut(p1); break;
    3417      116947 :     case 7: *s=subsr(-1,p1); *c=mpaut(p1); togglesign(*c); break;
    3418             :   }
    3419     1674608 :   gptr[0]=s; gptr[1]=c;
    3420     1674608 :   gerepilemanysp(av,tetpil,gptr,2);
    3421             : }
    3422             : 
    3423             : /* SINE and COSINE - 1 */
    3424             : void
    3425        4910 : mpsincosm1(GEN x, GEN *s, GEN *c)
    3426             : {
    3427             :   long mod8;
    3428             :   pari_sp av, tetpil;
    3429             :   GEN p1, *gptr[2];
    3430             : 
    3431        4910 :   if (!signe(x))
    3432             :   {
    3433           0 :     long e = expo(x);
    3434           0 :     *s = real_0_bit(e);
    3435           0 :     *c = real_0_bit(2*e-1);
    3436           0 :     return;
    3437             :   }
    3438        4910 :   av=avma; p1=mpcosm1(x,&mod8); tetpil=avma;
    3439        4910 :   switch(mod8)
    3440             :   {
    3441        4490 :     case 0: *c=rcopy(p1); *s=mpaut(p1); break;
    3442          42 :     case 1: *s=addsr(1,p1); *c=addrs(mpaut(p1),1); togglesign(*c); break;
    3443           0 :     case 2: *c=subsr(-2,p1); *s=mpaut(p1); togglesign(*s); break;
    3444           0 :     case 3: *s=subsr(-1,p1); *c=subrs(mpaut(p1),1); break;
    3445         287 :     case 4: *c=rcopy(p1); *s=mpaut(p1); togglesign(*s); break;
    3446          77 :     case 5: *s=addsr( 1,p1); *c=subrs(mpaut(p1),1); break;
    3447           7 :     case 6: *c=subsr(-2,p1); *s=mpaut(p1); break;
    3448           7 :     case 7: *s=subsr(-1,p1); *c=subsr(-1,mpaut(p1)); break;
    3449             :   }
    3450        4910 :   gptr[0]=s; gptr[1]=c;
    3451        4910 :   gerepilemanysp(av,tetpil,gptr,2);
    3452             : }
    3453             : 
    3454             : /* return exp(ix), x a t_REAL */
    3455             : GEN
    3456      132277 : expIr(GEN x)
    3457             : {
    3458      132277 :   pari_sp av = avma;
    3459      132277 :   GEN v = cgetg(3,t_COMPLEX);
    3460      132277 :   mpsincos(x, (GEN*)(v+2), (GEN*)(v+1));
    3461      132277 :   if (!signe(gel(v,2))) return gerepilecopy(av, gel(v,1));
    3462      131731 :   return v;
    3463             : }
    3464             : 
    3465             : /* return exp(ix)-1, x a t_REAL */
    3466             : static GEN
    3467        4910 : expm1_Ir(GEN x)
    3468             : {
    3469        4910 :   pari_sp av = avma;
    3470        4910 :   GEN v = cgetg(3,t_COMPLEX);
    3471        4910 :   mpsincosm1(x, (GEN*)(v+2), (GEN*)(v+1));
    3472        4910 :   if (!signe(gel(v,2))) return gerepilecopy(av, gel(v,1));
    3473        4910 :   return v;
    3474             : }
    3475             : 
    3476             : /* return exp(z)-1, z complex */
    3477             : GEN
    3478        5008 : cxexpm1(GEN z, long prec)
    3479             : {
    3480        5008 :   pari_sp av = avma;
    3481        5008 :   GEN X, Y, x = real_i(z), y = imag_i(z);
    3482        5008 :   long l = precision(z);
    3483        5008 :   if (l) prec = l;
    3484        5008 :   if (typ(x) != t_REAL) x = gtofp(x, prec);
    3485        5008 :   if (typ(y) != t_REAL) y = gtofp(y, prec);
    3486        5008 :   if (gequal0(y)) return mpexpm1(x);
    3487        4910 :   if (gequal0(x)) return expm1_Ir(y);
    3488        4826 :   X = mpexpm1(x); /* t_REAL */
    3489        4826 :   Y = expm1_Ir(y);
    3490             :   /* exp(x+iy) - 1 = (exp(x)-1)(exp(iy)-1) + exp(x)-1 + exp(iy)-1 */
    3491        4826 :   return gerepileupto(av, gadd(gadd(X,Y), gmul(X,Y)));
    3492             : }
    3493             : 
    3494             : void
    3495     1423350 : gsincos(GEN x, GEN *s, GEN *c, long prec)
    3496             : {
    3497             :   long i, j, ex, ex2, lx, ly, mi;
    3498             :   pari_sp av, tetpil;
    3499             :   GEN y, r, u, v, u1, v1, p1, p2, p3, p4, ps, pc;
    3500             :   GEN *gptr[4];
    3501             : 
    3502     1423350 :   switch(typ(x))
    3503             :   {
    3504             :     case t_INT: case t_FRAC:
    3505         959 :       *s = cgetr(prec);
    3506         959 :       *c = cgetr(prec); av = avma;
    3507         959 :       mpsincos(tofp_safe(x, prec), &ps, &pc);
    3508         959 :       affrr_fixlg(ps,*s);
    3509     1424309 :       affrr_fixlg(pc,*c); set_avma(av); return;
    3510             : 
    3511             :     case t_REAL:
    3512     1417876 :       mpsincos(x,s,c); return;
    3513             : 
    3514             :     case t_COMPLEX:
    3515        4116 :       i = precision(x); if (i) prec = i;
    3516        4116 :       ps = cgetc(prec); *s = ps;
    3517        4116 :       pc = cgetc(prec); *c = pc; av = avma;
    3518        4116 :       r = gexp(gel(x,2),prec);
    3519        4116 :       v1 = gmul2n(addrr(invr(r),r), -1); /* = cos(I*Im(x)) */
    3520        4116 :       u1 = subrr(r, v1); /* = I*sin(I*Im(x)) */
    3521        4116 :       gsincos(gel(x,1), &u,&v, prec);
    3522        4116 :       affrr_fixlg(mulrr(v1,u), gel(ps,1));
    3523        4116 :       affrr_fixlg(mulrr(u1,v), gel(ps,2));
    3524        4116 :       affrr_fixlg(mulrr(v1,v), gel(pc,1));
    3525        4116 :       affrr_fixlg(mulrr(u1,u), gel(pc,2)); togglesign(gel(pc,2));
    3526        4116 :       set_avma(av); return;
    3527             : 
    3528             :     case t_QUAD:
    3529           0 :       av = avma; gsincos(quadtofp(x, prec), s, c, prec);
    3530           0 :       gerepileall(av, 2, s, c); return;
    3531             : 
    3532             :     default:
    3533         399 :       av = avma; if (!(y = toser_i(x))) break;
    3534         399 :       if (gequal0(y)) { *s = gerepilecopy(av,y); *c = gaddsg(1,*s); return; }
    3535             : 
    3536         399 :       ex = valp(y); lx = lg(y); ex2 = 2*ex+2;
    3537         399 :       if (ex < 0) pari_err_DOMAIN("gsincos","valuation", "<", gen_0, x);
    3538         399 :       if (ex2 > lx)
    3539             :       {
    3540         105 :         *s = x == y? gcopy(y): gerepilecopy(av, y); av = avma;
    3541         105 :         *c = gerepileupto(av, gsubsg(1, gdivgs(gsqr(y),2)));
    3542         105 :         return;
    3543             :       }
    3544         294 :       if (!ex)
    3545             :       {
    3546          63 :         gsincos(serchop0(y),&u,&v,prec);
    3547          63 :         gsincos(gel(y,2),&u1,&v1,prec);
    3548          63 :         p1 = gmul(v1,v);
    3549          63 :         p2 = gmul(u1,u);
    3550          63 :         p3 = gmul(v1,u);
    3551          63 :         p4 = gmul(u1,v); tetpil = avma;
    3552          63 :         *c = gsub(p1,p2);
    3553          63 :         *s = gadd(p3,p4);
    3554          63 :         gptr[0]=s; gptr[1]=c;
    3555          63 :         gerepilemanysp(av,tetpil,gptr,2);
    3556          63 :         return;
    3557             :       }
    3558             : 
    3559         231 :       ly = lx+2*ex;
    3560         231 :       mi = lx-1; while (mi>=3 && isrationalzero(gel(y,mi))) mi--;
    3561         231 :       mi += ex-2;
    3562         231 :       pc = cgetg(ly,t_SER); *c = pc;
    3563         231 :       ps = cgetg(lx,t_SER); *s = ps;
    3564         231 :       pc[1] = evalsigne(1) | _evalvalp(0) | evalvarn(varn(y));
    3565         231 :       gel(pc,2) = gen_1; ps[1] = y[1];
    3566         231 :       for (i=2; i<ex+2; i++) gel(ps,i) = gcopy(gel(y,i));
    3567         231 :       for (i=3; i< ex2; i++) gel(pc,i) = gen_0;
    3568        3129 :       for (i=ex2; i<ly; i++)
    3569             :       {
    3570        2898 :         long ii = i-ex;
    3571        2898 :         av = avma; p1 = gen_0;
    3572        6678 :         for (j=ex; j<=minss(ii-2,mi); j++)
    3573        3780 :           p1 = gadd(p1, gmulgs(gmul(gel(y,j-ex+2),gel(ps,ii-j)),j));
    3574        2898 :         gel(pc,i) = gerepileupto(av, gdivgs(p1,2-i));
    3575        2898 :         if (ii < lx)
    3576             :         {
    3577        2660 :           av = avma; p1 = gen_0;
    3578        5726 :           for (j=ex; j<=minss(i-ex2,mi); j++)
    3579        3066 :             p1 = gadd(p1,gmulgs(gmul(gel(y,j-ex+2),gel(pc,i-j)),j));
    3580        2660 :           p1 = gdivgs(p1,i-2);
    3581        2660 :           gel(ps,ii) = gerepileupto(av, gadd(p1,gel(y,ii)));
    3582             :         }
    3583             :       }
    3584         231 :       return;
    3585             :   }
    3586           0 :   pari_err_TYPE("gsincos",x);
    3587             : }
    3588             : 
    3589             : /********************************************************************/
    3590             : /**                                                                **/
    3591             : /**                              SINC                              **/
    3592             : /**                                                                **/
    3593             : /********************************************************************/
    3594             : static GEN
    3595      107450 : mpsinc(GEN x)
    3596             : {
    3597      107450 :   pari_sp av = avma;
    3598             :   GEN s, c;
    3599             : 
    3600      107450 :   if (!signe(x)) {
    3601           0 :     long l = nbits2prec(-expo(x));
    3602           0 :     if (l < LOWDEFAULTPREC) l = LOWDEFAULTPREC;
    3603           0 :     return real_1(l);
    3604             :   }
    3605             : 
    3606      107450 :   mpsincos(x,&s,&c);
    3607      107450 :   return gerepileuptoleaf(av, divrr(s,x));
    3608             : }
    3609             : 
    3610             : GEN
    3611      107562 : gsinc(GEN x, long prec)
    3612             : {
    3613             :   pari_sp av;
    3614             :   GEN r, u, v, y, u1, v1;
    3615             :   long i;
    3616             : 
    3617      107562 :   switch(typ(x))
    3618             :   {
    3619      107429 :     case t_REAL: return mpsinc(x);
    3620             :     case t_COMPLEX:
    3621          49 :       if (isintzero(gel(x,1)))
    3622             :       {
    3623          28 :         av = avma; x = gel(x,2);
    3624          28 :         if (gequal0(x)) return gcosh(x,prec);
    3625          14 :         return gerepileuptoleaf(av,gdiv(gsinh(x,prec),x));
    3626             :       }
    3627          21 :       i = precision(x); if (i) prec = i;
    3628          21 :       y = cgetc(prec); av = avma;
    3629          21 :       r = gexp(gel(x,2),prec);
    3630          21 :       v1 = gmul2n(addrr(invr(r),r), -1); /* = cos(I*Im(x)) */
    3631          21 :       u1 = subrr(r, v1); /* = I*sin(I*Im(x)) */
    3632          21 :       gsincos(gel(x,1),&u,&v,prec);
    3633          21 :       affc_fixlg(gdiv(mkcomplex(gmul(v1,u), gmul(u1,v)), x), y);
    3634          21 :       set_avma(av); return y;
    3635             : 
    3636             :     case t_INT:
    3637          14 :       if (!signe(x)) return real_1(prec); /*fall through*/
    3638             :     case t_FRAC:
    3639          21 :       y = cgetr(prec); av = avma;
    3640          21 :       affrr_fixlg(mpsinc(tofp_safe(x,prec)), y); set_avma(av); return y;
    3641             : 
    3642             :     case t_PADIC:
    3643          21 :       if (gequal0(x)) return cvtop(gen_1, gel(x,2), valp(x));
    3644          14 :       av = avma; y = sin_p(x);
    3645          14 :       if (!y) pari_err_DOMAIN("gsinc(t_PADIC)","argument","",gen_0,x);
    3646           7 :       return gerepileupto(av,gdiv(y,x));
    3647             : 
    3648             :     default:
    3649             :     {
    3650             :       long ex;
    3651          35 :       av = avma; if (!(y = toser_i(x))) break;
    3652          35 :       if (gequal0(y)) return gerepileupto(av, gaddsg(1,y));
    3653          35 :       ex = valp(y);
    3654          35 :       if (ex < 0) pari_err_DOMAIN("sinc","valuation", "<", gen_0, x);
    3655          28 :       if (ex)
    3656             :       {
    3657          28 :         gsincos(y,&u,&v,prec);
    3658          28 :         y = gerepileupto(av, gdiv(u,y));
    3659          28 :         if (lg(y) > 2) gel(y,2) = gen_1;
    3660          28 :         return y;
    3661             :       }
    3662             :       else
    3663             :       {
    3664           0 :         GEN z0, y0 = gel(y,2), y1 = serchop0(y), y10 = y1;
    3665           0 :         if (!gequal1(y0)) y10 = gdiv(y10, y0);
    3666           0 :         gsincos(y1,&u,&v,prec);
    3667           0 :         z0 = gdiv(gcos(y0,prec), y0);
    3668           0 :         y = gaddsg(1, y10);
    3669           0 :         u = gadd(gmul(gsinc(y0, prec),v), gmul(z0, u));
    3670           0 :         return gerepileupto(av,gdiv(u,y));
    3671             :       }
    3672             :     }
    3673             :   }
    3674           0 :   return trans_eval("sinc",gsinc,x,prec);
    3675             : }
    3676             : 
    3677             : /********************************************************************/
    3678             : /**                                                                **/
    3679             : /**                     TANGENT and COTANGENT                      **/
    3680             : /**                                                                **/
    3681             : /********************************************************************/
    3682             : static GEN
    3683          35 : mptan(GEN x)
    3684             : {
    3685          35 :   pari_sp av = avma;
    3686             :   GEN s, c;
    3687             : 
    3688          35 :   mpsincos(x,&s,&c);
    3689          35 :   if (!signe(c))
    3690           0 :     pari_err_DOMAIN("tan", "argument", "=", strtoGENstr("Pi/2 + kPi"),x);
    3691          35 :   return gerepileuptoleaf(av, divrr(s,c));
    3692             : }
    3693             : 
    3694             : GEN
    3695         105 : gtan(GEN x, long prec)
    3696             : {
    3697             :   pari_sp av;
    3698             :   GEN y, s, c;
    3699             : 
    3700         105 :   switch(typ(x))
    3701             :   {
    3702          28 :     case t_REAL: return mptan(x);
    3703             : 
    3704             :     case t_COMPLEX: {
    3705          28 :       if (isintzero(gel(x,1))) retmkcomplex(gen_0,gtanh(gel(x,2),prec));
    3706          14 :       av = avma; y = mulcxmI(gtanh(mulcxI(x), prec)); /* tan x = -I th(I x) */
    3707          14 :       gel(y,1) = gcopy(gel(y,1));
    3708          14 :       return gerepileupto(av, y);
    3709             :     }
    3710             :     case t_INT: case t_FRAC:
    3711           7 :       y = cgetr(prec); av = avma;
    3712           7 :       affrr_fixlg(mptan(tofp_safe(x,prec)), y); set_avma(av); return y;
    3713             : 
    3714             :     case t_PADIC:
    3715          14 :       av = avma;
    3716          14 :       return gerepileupto(av, gdiv(gsin(x,prec), gcos(x,prec)));
    3717             : 
    3718             :     default:
    3719          28 :       av = avma; if (!(y = toser_i(x))) break;
    3720          21 :       if (gequal0(y)) return gerepilecopy(av, y);
    3721          21 :       if (valp(y) < 0)
    3722           7 :         pari_err_DOMAIN("tan","valuation", "<", gen_0, x);
    3723          14 :       gsincos(y,&s,&c,prec);
    3724          14 :       return gerepileupto(av, gdiv(s,c));
    3725             :   }
    3726           7 :   return trans_eval("tan",gtan,x,prec);
    3727             : }
    3728             : 
    3729             : static GEN
    3730          63 : mpcotan(GEN x)
    3731             : {
    3732          63 :   pari_sp av=avma, tetpil;
    3733             :   GEN s,c;
    3734             : 
    3735          63 :   mpsincos(x,&s,&c); tetpil=avma;
    3736          63 :   return gerepile(av,tetpil,divrr(c,s));
    3737             : }
    3738             : 
    3739             : GEN
    3740        4151 : gcotan(GEN x, long prec)
    3741             : {
    3742             :   pari_sp av;
    3743             :   GEN y, s, c;
    3744             : 
    3745        4151 :   switch(typ(x))
    3746             :   {
    3747             :     case t_REAL:
    3748          56 :       return mpcotan(x);
    3749             : 
    3750             :     case t_COMPLEX:
    3751        4004 :       if (isintzero(gel(x,1))) {
    3752          21 :         GEN z = cgetg(3, t_COMPLEX);
    3753          21 :         gel(z,1) = gen_0;
    3754          21 :         av = avma;
    3755          21 :         gel(z,2) = gerepileupto(av, gneg(ginv(gtanh(gel(x,2),prec))));
    3756          21 :         return z;
    3757             :       }
    3758        3983 :       av = avma;
    3759        3983 :       gsincos(x,&s,&c,prec);
    3760        3983 :       return gerepileupto(av, gdiv(c,s));
    3761             : 
    3762             :     case t_INT: case t_FRAC:
    3763           7 :       y = cgetr(prec); av = avma;
    3764           7 :       affrr_fixlg(mpcotan(tofp_safe(x,prec)), y); set_avma(av); return y;
    3765             : 
    3766             :     case t_PADIC:
    3767          14 :       av = avma;
    3768          14 :       return gerepileupto(av, gdiv(gcos(x,prec), gsin(x,prec)));
    3769             : 
    3770             :     default:
    3771          70 :       av = avma; if (!(y = toser_i(x))) break;
    3772          63 :       if (gequal0(y)) pari_err_DOMAIN("cotan", "argument", "=", gen_0, y);
    3773          63 :       if (valp(y) < 0) pari_err_DOMAIN("cotan","valuation", "<", gen_0, x);
    3774          56 :       gsincos(y,&s,&c,prec);
    3775          56 :       return gerepileupto(av, gdiv(c,s));
    3776             :   }
    3777           7 :   return trans_eval("cotan",gcotan,x,prec);
    3778             : }

Generated by: LCOV version 1.13