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 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.10.0 lcov report (development 20277-2bd9113) Lines: 1938 1988 97.5 %
Date: 2017-02-19 05:49:50 Functions: 141 143 98.6 %
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       64128 : pari_init_floats(void)
      31             : {
      32       64128 :   gcatalan = geuler = gpi = bernzone = glog2 = NULL;
      33       64128 : }
      34             : 
      35             : void
      36       64448 : pari_close_floats(void)
      37             : {
      38       64448 :   if (gcatalan) gunclone(gcatalan);
      39       64568 :   if (geuler) gunclone(geuler);
      40       64568 :   if (gpi) gunclone(gpi);
      41       64568 :   if (bernzone) gunclone(bernzone);
      42       64568 :   if (glog2) gunclone(glog2);
      43       64568 : }
      44             : 
      45             : /********************************************************************/
      46             : /**                   GENERIC BINARY SPLITTING                     **/
      47             : /**                    (Haible, Papanikolaou)                      **/
      48             : /********************************************************************/
      49             : void
      50        6505 : abpq_init(struct abpq *A, long n)
      51             : {
      52        6505 :   A->a = (GEN*)new_chunk(n+1);
      53        6505 :   A->b = (GEN*)new_chunk(n+1);
      54        6505 :   A->p = (GEN*)new_chunk(n+1);
      55        6505 :   A->q = (GEN*)new_chunk(n+1);
      56        6505 : }
      57             : static GEN
      58      567479 : mulii3(GEN a, GEN b, GEN c) { return mulii(mulii(a,b),c); }
      59             : static GEN
      60      168302 : 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      168302 : T2(struct abpq *A, long n1, GEN P)
      65             : {
      66      168302 :   GEN u1 = mulii4(A->a[n1], A->p[n1], A->b[n1+1], A->q[n1+1]);
      67      168302 :   GEN u2 = mulii3(A->b[n1],A->a[n1+1], P);
      68      168302 :   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      330171 : 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      330171 :   switch(n2 - n1)
      80             :   {
      81             :     GEN b, p, q;
      82             :     case 1:
      83          36 :       r->P = A->p[n1];
      84          36 :       r->Q = A->q[n1];
      85          36 :       r->B = A->b[n1];
      86          36 :       r->T = mulii(A->a[n1], A->p[n1]);
      87      168374 :       return;
      88             :     case 2:
      89       92791 :       r->P = mulii(A->p[n1], A->p[n1+1]);
      90       92791 :       r->Q = mulii(A->q[n1], A->q[n1+1]);
      91       92791 :       r->B = mulii(A->b[n1], A->b[n1+1]);
      92       92791 :       av = avma;
      93       92791 :       r->T = gerepileuptoint(av, T2(A, n1, r->P));
      94       92791 :       return;
      95             : 
      96             :     case 3:
      97       75511 :       p = mulii(A->p[n1+1], A->p[n1+2]);
      98       75511 :       q = mulii(A->q[n1+1], A->q[n1+2]);
      99       75511 :       b = mulii(A->b[n1+1], A->b[n1+2]);
     100       75511 :       r->P = mulii(A->p[n1], p);
     101       75511 :       r->Q = mulii(A->q[n1], q);
     102       75511 :       r->B = mulii(A->b[n1], b);
     103       75511 :       av = avma;
     104       75511 :       u1 = mulii3(b, q, A->a[n1]);
     105       75511 :       u2 = mulii(A->b[n1], T2(A, n1+1, p));
     106       75511 :       r->T = gerepileuptoint(av, mulii(A->p[n1], addii(u1, u2)));
     107       75511 :       return;
     108             :   }
     109             : 
     110      161833 :   av = avma;
     111      161833 :   n = (n1 + n2) >> 1;
     112      161833 :   abpq_sum(&L, n1, n, A);
     113      161833 :   abpq_sum(&R, n, n2, A);
     114             : 
     115      161833 :   r->P = mulii(L.P, R.P);
     116      161833 :   r->Q = mulii(L.Q, R.Q);
     117      161833 :   r->B = mulii(L.B, R.B);
     118      161833 :   u1 = mulii3(R.B,R.Q,L.T);
     119      161833 :   u2 = mulii3(L.B,L.P,R.T);
     120      161833 :   r->T = addii(u1,u2);
     121      161833 :   avma = av;
     122      161833 :   r->P = icopy(r->P);
     123      161833 :   r->Q = icopy(r->Q);
     124      161833 :   r->B = icopy(r->B);
     125      161833 :   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        4642 : swap_clone(GEN *old, GEN c)
     136             : {
     137        4642 :   GEN tmp = *old;
     138        4642 :   *old = c;
     139        4642 :   if (tmp) gunclone(tmp);
     140        4642 : }
     141             : 
     142             : /*                         ----
     143             :  *  53360 (640320)^(1/2)   \    (6n)! (545140134 n + 13591409)
     144             :  *  -------------------- = /    ------------------------------
     145             :  *        Pi               ----   (n!)^3 (3n)! (-640320)^(3n)
     146             :  *                         n>=0
     147             :  *
     148             :  * Ramanujan's formula + binary splitting */
     149             : static GEN
     150        2057 : pi_ramanujan(long prec)
     151             : {
     152        2057 :   const ulong B = 545140134, A = 13591409, C = 640320;
     153        2057 :   const double alpha2 = 47.11041314; /* 3log(C/12) / log(2) */
     154             :   long n, nmax, prec2;
     155             :   struct abpq_res R;
     156             :   struct abpq S;
     157             :   GEN D, u;
     158             : 
     159        2057 :   nmax = (long)(1 + prec2nbits(prec)/alpha2);
     160             : #ifdef LONG_IS_64BIT
     161        1748 :   D = utoipos(10939058860032000UL); /* C^3/24 */
     162             : #else
     163         309 :   D = uutoi(2546948UL,495419392UL);
     164             : #endif
     165        2057 :   abpq_init(&S, nmax);
     166        2057 :   S.a[0] = utoipos(A);
     167        2057 :   S.b[0] = S.p[0] = S.q[0] = gen_1;
     168       35756 :   for (n = 1; n <= nmax; n++)
     169             :   {
     170       33699 :     S.a[n] = addiu(muluu(B, n), A);
     171       33699 :     S.b[n] = gen_1;
     172       33699 :     S.p[n] = mulis(muluu(6*n-5, 2*n-1), 1-6*n);
     173       33699 :     S.q[n] = mulii(sqru(n), muliu(D,n));
     174             :   }
     175        2057 :   abpq_sum(&R, 0, nmax, &S); prec2 = prec+EXTRAPRECWORD;
     176        2057 :   u = itor(muliu(R.Q,C/12), prec2);
     177        2057 :   return rtor(mulrr(divri(u, R.T), sqrtr_abs(utor(C,prec2))), prec);
     178             : }
     179             : 
     180             : #if 0 /* Much slower than binary splitting at least up to prec = 10^8 */
     181             : /* Gauss - Brent-Salamin AGM iteration */
     182             : static GEN
     183             : pi_brent_salamin(long prec)
     184             : {
     185             :   GEN A, B, C;
     186             :   pari_sp av2;
     187             :   long i, G;
     188             : 
     189             :   G = - prec2nbits(prec);
     190             :   incrprec(prec);
     191             : 
     192             :   A = real2n(-1, prec);
     193             :   B = sqrtr_abs(A); /* = 1/sqrt(2) */
     194             :   setexpo(A, 0);
     195             :   C = real2n(-2, prec); av2 = avma;
     196             :   for (i = 0;; i++)
     197             :   {
     198             :     GEN y, a, b, B_A = subrr(B, A);
     199             :     pari_sp av3 = avma;
     200             :     if (expo(B_A) < G) break;
     201             :     a = addrr(A,B); shiftr_inplace(a, -1);
     202             :     b = mulrr(A,B);
     203             :     affrr(a, A);
     204             :     affrr(sqrtr_abs(b), B); avma = av3;
     205             :     y = sqrr(B_A); shiftr_inplace(y, i - 2);
     206             :     affrr(subrr(C, y), C); avma = av2;
     207             :   }
     208             :   shiftr_inplace(C, 2);
     209             :   return divrr(sqrr(addrr(A,B)), C);
     210             : }
     211             : GEN
     212             : constpi(long prec)
     213             : {
     214             :   pari_sp av;
     215             :   GEN tmp;
     216             :   if (gpi && realprec(gpi) >= prec) return gpi;
     217             : 
     218             :   tmp = cgetr_block(prec);
     219             :   av = avma;
     220             :   affrr(pi_brent_salamin(prec), tmp);
     221             :   swap_clone(&gpi, tmp);
     222             :   avma = av;  return gpi;
     223             : }
     224             : #endif
     225             : 
     226             : GEN
     227    10506246 : constpi(long prec)
     228             : {
     229             :   pari_sp av;
     230             :   GEN tmp;
     231    10506246 :   if (gpi && realprec(gpi) >= prec) return gpi;
     232             : 
     233        2057 :   av = avma;
     234        2057 :   tmp = gclone(pi_ramanujan(prec));
     235        2057 :   swap_clone(&gpi,tmp);
     236        2057 :   avma = av; return gpi;
     237             : }
     238             : 
     239             : GEN
     240    10505555 : mppi(long prec) { return rtor(constpi(prec), prec); }
     241             : 
     242             : /* Pi * 2^n */
     243             : GEN
     244     4354218 : Pi2n(long n, long prec)
     245             : {
     246     4354218 :   GEN x = mppi(prec); shiftr_inplace(x, n);
     247     4354218 :   return x;
     248             : }
     249             : 
     250             : /* I * Pi * 2^n */
     251             : GEN
     252        6166 : PiI2n(long n, long prec) { retmkcomplex(gen_0, Pi2n(n, prec)); }
     253             : 
     254             : /* 2I * Pi */
     255             : GEN
     256        4354 : PiI2(long prec) { return PiI2n(1, prec); }
     257             : 
     258             : /********************************************************************/
     259             : /**                                                                **/
     260             : /**                       EULER CONSTANT                           **/
     261             : /**                                                                **/
     262             : /********************************************************************/
     263             : 
     264             : GEN
     265       19792 : consteuler(long prec)
     266             : {
     267             :   GEN u,v,a,b,tmpeuler;
     268             :   long l, n1, n, k, x;
     269             :   pari_sp av1, av2;
     270             : 
     271       19792 :   if (geuler && realprec(geuler) >= prec) return geuler;
     272             : 
     273         354 :   av1 = avma; tmpeuler = cgetr_block(prec);
     274             : 
     275         354 :   incrprec(prec);
     276             : 
     277         354 :   l = prec+EXTRAPRECWORD; x = (long) (1 + prec2nbits_mul(l, LOG2/4));
     278         354 :   a = utor(x,l); u=logr_abs(a); setsigne(u,-1); affrr(u,a);
     279         354 :   b = real_1(l);
     280         354 :   v = real_1(l);
     281         354 :   n = (long)(1+3.591*x); /* z=3.591: z*[ ln(z)-1 ]=1 */
     282         354 :   n1 = minss(n, SQRTVERYBIGINT);
     283         354 :   if (x < SQRTVERYBIGINT)
     284             :   {
     285         354 :     ulong xx = x*x;
     286         354 :     av2 = avma;
     287      126187 :     for (k=1; k<n1; k++)
     288             :     {
     289      125833 :       affrr(divru(mulur(xx,b),k*k), b);
     290      125833 :       affrr(divru(addrr(divru(mulur(xx,a),k),b),k), a);
     291      125833 :       affrr(addrr(u,a), u);
     292      125833 :       affrr(addrr(v,b), v); avma = av2;
     293             :     }
     294         708 :     for (   ; k<=n; k++)
     295             :     {
     296         354 :       affrr(divru(divru(mulur(xx,b),k),k), b);
     297         354 :       affrr(divru(addrr(divru(mulur(xx,a),k),b),k), a);
     298         354 :       affrr(addrr(u,a), u);
     299         354 :       affrr(addrr(v,b), v); avma = av2;
     300             :     }
     301             :   }
     302             :   else
     303             :   {
     304           0 :     GEN xx = sqru(x);
     305           0 :     av2 = avma;
     306           0 :     for (k=1; k<n1; k++)
     307             :     {
     308           0 :       affrr(divru(mulir(xx,b),k*k), b);
     309           0 :       affrr(divru(addrr(divru(mulir(xx,a),k),b),k), a);
     310           0 :       affrr(addrr(u,a), u);
     311           0 :       affrr(addrr(v,b), v); avma = av2;
     312             :     }
     313           0 :     for (   ; k<=n; k++)
     314             :     {
     315           0 :       affrr(divru(divru(mulir(xx,b),k),k), b);
     316           0 :       affrr(divru(addrr(divru(mulir(xx,a),k),b),k), a);
     317           0 :       affrr(addrr(u,a), u);
     318           0 :       affrr(addrr(v,b), v); avma = av2;
     319             :     }
     320             :   }
     321         354 :   divrrz(u,v,tmpeuler);
     322         354 :   swap_clone(&geuler,tmpeuler);
     323         354 :   avma = av1; return geuler;
     324             : }
     325             : 
     326             : GEN
     327       19792 : mpeuler(long prec) { return rtor(consteuler(prec), prec); }
     328             : 
     329             : /********************************************************************/
     330             : /**                                                                **/
     331             : /**                       CATALAN CONSTANT                         **/
     332             : /**                                                                **/
     333             : /********************************************************************/
     334             : /* 8G = 3\sum_{n>=0} 1/(binomial(2n,n)(2n+1)^2) + Pi log(2+sqrt(3)) */
     335             : static GEN
     336          14 : catalan(long prec)
     337             : {
     338          14 :   long i, nmax = prec2nbits(prec) >> 1;
     339             :   struct abpq_res R;
     340             :   struct abpq A;
     341             :   GEN u, v;
     342          14 :   abpq_init(&A, nmax);
     343          14 :   A.a[0] = A.b[0] = A.p[0] = A.q[0] = gen_1;
     344        6510 :   for (i = 1; i <= nmax; i++)
     345             :   {
     346        6496 :     A.a[i] = gen_1;
     347        6496 :     A.b[i] = utoipos((i<<1)+1);
     348        6496 :     A.p[i] = utoipos(i);
     349        6496 :     A.q[i] = utoipos((i<<2)+2);
     350             :   }
     351          14 :   abpq_sum(&R, 0, nmax, &A);
     352          14 :   u = mulur(3, rdivii(R.T, mulii(R.B,R.Q),prec));
     353          14 :   v = mulrr(mppi(prec), logr_abs(addrs(sqrtr_abs(utor(3,prec)), 2)));
     354          14 :   u = addrr(u, v); shiftr_inplace(u, -3);
     355          14 :   return u;
     356             : }
     357             : 
     358             : GEN
     359          14 : constcatalan(long prec)
     360             : {
     361          14 :   pari_sp av = avma;
     362             :   GEN tmp;
     363          14 :   if (gcatalan && realprec(gcatalan) >= prec) return gcatalan;
     364          14 :   tmp = gclone(catalan(prec));
     365          14 :   swap_clone(&gcatalan,tmp);
     366          14 :   avma = av; return gcatalan;
     367             : }
     368             : 
     369             : GEN
     370          14 : mpcatalan(long prec) { return rtor(constcatalan(prec), prec); }
     371             : 
     372             : /********************************************************************/
     373             : /**                                                                **/
     374             : /**          TYPE CONVERSION FOR TRANSCENDENTAL FUNCTIONS          **/
     375             : /**                                                                **/
     376             : /********************************************************************/
     377             : static GEN
     378       41739 : transvec(GEN (*f)(GEN,long), GEN x, long prec)
     379             : {
     380             :   long i, l;
     381       41739 :   GEN y = cgetg_copy(x, &l);
     382       41739 :   for (i=1; i<l; i++) gel(y,i) = f(gel(x,i),prec);
     383       41732 :   return y;
     384             : }
     385             : 
     386             : GEN
     387      430154 : trans_eval(const char *fun, GEN (*f)(GEN,long), GEN x, long prec)
     388             : {
     389      430154 :   pari_sp av = avma;
     390      430154 :   if (prec < 3) pari_err_BUG("trans_eval [prec < 3]");
     391      430154 :   switch(typ(x))
     392             :   {
     393      373006 :     case t_INT:    x = f(itor(x,prec),prec); break;
     394       15353 :     case t_FRAC:   x = f(fractor(x, prec),prec); break;
     395           7 :     case t_QUAD:   x = f(quadtofp(x,prec),prec); break;
     396           7 :     case t_POLMOD: x = transvec(f, polmod_to_embed(x,prec), prec); break;
     397             :     case t_VEC:
     398             :     case t_COL:
     399       41732 :     case t_MAT: return transvec(f, x, prec);
     400          49 :     default: pari_err_TYPE(fun,x); return NULL;
     401             :   }
     402      388352 :   return gerepileupto(av, x);
     403             : }
     404             : 
     405             : /*******************************************************************/
     406             : /*                                                                 */
     407             : /*                            POWERING                             */
     408             : /*                                                                 */
     409             : /*******************************************************************/
     410             : /* x a t_REAL 0, return exp(x) */
     411             : static GEN
     412      235505 : mpexp0(GEN x)
     413             : {
     414      235505 :   long e = expo(x);
     415      235505 :   return e >= 0? real_0_bit(e): real_1_bit(-e);
     416             : }
     417             : static GEN
     418        1701 : powr0(GEN x)
     419        1701 : { return signe(x)? real_1(realprec(x)): mpexp0(x); }
     420             : 
     421             : /* x t_POL or t_SER, return scalarpol(RgX_get_1(x)) */
     422             : static GEN
     423      115619 : scalarpol_get_1(GEN x)
     424             : {
     425      115619 :   GEN y = cgetg(3,t_POL);
     426      115619 :   y[1] = evalvarn(varn(x)) | evalsigne(1);
     427      115619 :   gel(y,2) = RgX_get_1(x); return y;
     428             : }
     429             : /* to be called by the generic function gpowgs(x,s) when s = 0 */
     430             : static GEN
     431     2846472 : gpowg0(GEN x)
     432             : {
     433             :   long lx, i;
     434             :   GEN y;
     435             : 
     436     2846472 :   switch(typ(x))
     437             :   {
     438             :     case t_INT: case t_REAL: case t_FRAC: case t_PADIC:
     439      513403 :       return gen_1;
     440             : 
     441           7 :     case t_QUAD: x++; /*fall through*/
     442             :     case t_COMPLEX: {
     443       13731 :       pari_sp av = avma;
     444       13731 :       GEN a = gpowg0(gel(x,1));
     445       13731 :       GEN b = gpowg0(gel(x,2));
     446       13731 :       if (a == gen_1) return b;
     447          14 :       if (b == gen_1) return a;
     448           7 :       return gerepileupto(av, gmul(a,b));
     449             :     }
     450             :     case t_INTMOD:
     451         105 :       y = cgetg(3,t_INTMOD);
     452         105 :       gel(y,1) = icopy(gel(x,1));
     453         105 :       gel(y,2) = gen_1; return y;
     454             : 
     455        4599 :     case t_FFELT: return FF_1(x);
     456             : 
     457             :     case t_POLMOD:
     458          42 :       y = cgetg(3,t_POLMOD);
     459          42 :       gel(y,1) = gcopy(gel(x,1));
     460          42 :       gel(y,2) = scalarpol_get_1(gel(x,1)); return y;
     461             : 
     462             :     case t_RFRAC:
     463           7 :       return scalarpol_get_1(gel(x,2));
     464             :     case t_POL: case t_SER:
     465      115570 :       return scalarpol_get_1(x);
     466             : 
     467             :     case t_MAT:
     468          77 :       lx=lg(x); if (lx==1) return cgetg(1,t_MAT);
     469          70 :       if (lx != lgcols(x)) pari_err_DIM("gpow");
     470          70 :       y = matid(lx-1);
     471          70 :       for (i=1; i<lx; i++) gcoeff(y,i,i) = gpowg0(gcoeff(x,i,i));
     472          70 :       return y;
     473          14 :     case t_QFR: return qfr_1(x);
     474     2198910 :     case t_QFI: return qfi_1(x);
     475           7 :     case t_VECSMALL: return identity_perm(lg(x) - 1);
     476             :   }
     477           7 :   pari_err_TYPE("gpow",x);
     478             :   return NULL; /* LCOV_EXCL_LINE */
     479             : }
     480             : 
     481             : static GEN
     482    35247900 : _sqr(void *data /* ignored */, GEN x) { (void)data; return gsqr(x); }
     483             : static GEN
     484    14202485 : _mul(void *data /* ignored */, GEN x, GEN y) { (void)data; return gmul(x,y); }
     485             : static GEN
     486       45846 : _one(void *x) { return gpowg0((GEN) x); }
     487             : static GEN
     488     4586113 : _sqri(void *data /* ignored */, GEN x) { (void)data; return sqri(x); }
     489             : static GEN
     490     3069990 : _muli(void *data /* ignored */, GEN x, GEN y) { (void)data; return mulii(x,y); }
     491             : static GEN
     492     7108027 : _sqrr(void *data /* ignored */, GEN x) { (void)data; return sqrr(x); }
     493             : static GEN
     494     1290786 : _mulr(void *data /* ignored */, GEN x, GEN y) { (void)data; return mulrr(x,y); }
     495             : static GEN
     496       24011 : _oner(void *data /* prec */) { return real_1( *(long*) data); }
     497             : 
     498             : /* INTEGER POWERING (a^n for integer a != 0 and integer n > 0)
     499             :  *
     500             :  * Use left shift binary algorithm (RS is wasteful: multiplies big numbers,
     501             :  * with LS one of them is the base, hence small). Sign of result is set
     502             :  * to s (= 1,-1). Makes life easier for caller, which otherwise might do a
     503             :  * setsigne(gen_1 / gen_m1) */
     504             : static GEN
     505    28258781 : powiu_sign(GEN a, ulong N, long s)
     506             : {
     507             :   pari_sp av;
     508             :   GEN y;
     509             : 
     510    28258781 :   if (lgefint(a) == 3)
     511             :   { /* easy if |a| < 3 */
     512    24477142 :     ulong q = a[2];
     513    24477142 :     if (q == 1) return (s>0)? gen_1: gen_m1;
     514    15990584 :     if (q == 2) { a = int2u(N); setsigne(a,s); return a; }
     515    10047011 :     q = upowuu(q, N);
     516    10047012 :     if (q) return s>0? utoipos(q): utoineg(q);
     517             :   }
     518     4325462 :   if (N <= 2) {
     519     2745981 :     if (N == 2) return sqri(a);
     520        3616 :     a = icopy(a); setsigne(a,s); return a;
     521             :   }
     522     1579481 :   av = avma;
     523     1579481 :   y = gen_powu_i(a, N, NULL, &_sqri, &_muli);
     524     1579481 :   setsigne(y,s); return gerepileuptoint(av, y);
     525             : }
     526             : /* a^n */
     527             : GEN
     528    28065975 : powiu(GEN a, ulong n)
     529             : {
     530             :   long s;
     531    28065975 :   if (!n) return gen_1;
     532    27715806 :   s = signe(a);
     533    27715806 :   if (!s) return gen_0;
     534    27703835 :   return powiu_sign(a, n, (s < 0 && odd(n))? -1: 1);
     535             : }
     536             : GEN
     537    15596980 : powis(GEN a, long n)
     538             : {
     539             :   long s;
     540             :   GEN t, y;
     541    15596980 :   if (n >= 0) return powiu(a, n);
     542      106838 :   s = signe(a);
     543      106838 :   if (!s) pari_err_INV("powis",gen_0);
     544      106838 :   t = (s < 0 && odd(n))? gen_m1: gen_1;
     545      106838 :   if (is_pm1(a)) return t;
     546             :   /* n < 0, |a| > 1 */
     547      106684 :   y = cgetg(3,t_FRAC);
     548      106684 :   gel(y,1) = t;
     549      106684 :   gel(y,2) = powiu_sign(a, -n, 1); /* force denominator > 0 */
     550      106684 :   return y;
     551             : }
     552             : GEN
     553    10098434 : powuu(ulong p, ulong N)
     554             : {
     555    10098434 :   pari_sp av = avma;
     556    10098434 :   long P[] = {evaltyp(t_INT)|_evallg(3), evalsigne(1)|evallgefint(3),0};
     557             :   ulong pN;
     558             :   GEN y;
     559    10098434 :   if (N <= 2)
     560             :   {
     561     6950513 :     if (N == 2) return sqru(p);
     562     4835307 :     if (N == 1) return utoipos(p);
     563     2491866 :     return gen_1;
     564             :   }
     565     3147921 :   if (!p) return gen_0;
     566     3147921 :   pN = upowuu(p, N);
     567     3147921 :   if (pN) return utoipos(pN);
     568      738709 :   if (p == 2) return int2u(N);
     569      737994 :   P[2] = p; av = avma;
     570      737994 :   y = gen_powu_i(P, N, NULL, &_sqri, &_muli);
     571      737994 :   return gerepileuptoint(av, y);
     572             : }
     573             : 
     574             : /* return 0 if overflow */
     575             : static ulong
     576     2716996 : usqru(ulong p) { return p & HIGHMASK? 0: p*p; }
     577             : ulong
     578    18256218 : upowuu(ulong p, ulong k)
     579             : {
     580             : #ifdef LONG_IS_64BIT
     581    15631390 :   const ulong CUTOFF3 = 2642245;
     582    15631390 :   const ulong CUTOFF4 = 65535;
     583    15631390 :   const ulong CUTOFF5 = 7131;
     584    15631390 :   const ulong CUTOFF6 = 1625;
     585    15631390 :   const ulong CUTOFF7 = 565;
     586    15631390 :   const ulong CUTOFF8 = 255;
     587    15631390 :   const ulong CUTOFF9 = 138;
     588    15631390 :   const ulong CUTOFF10 = 84;
     589    15631390 :   const ulong CUTOFF11 = 56;
     590    15631390 :   const ulong CUTOFF12 = 40;
     591    15631390 :   const ulong CUTOFF13 = 30;
     592    15631390 :   const ulong CUTOFF14 = 23;
     593    15631390 :   const ulong CUTOFF15 = 19;
     594    15631390 :   const ulong CUTOFF16 = 15;
     595    15631390 :   const ulong CUTOFF17 = 13;
     596    15631390 :   const ulong CUTOFF18 = 11;
     597    15631390 :   const ulong CUTOFF19 = 10;
     598    15631390 :   const ulong CUTOFF20 =  9;
     599             : #else
     600     2624828 :   const ulong CUTOFF3 = 1625;
     601     2624828 :   const ulong CUTOFF4 =  255;
     602     2624828 :   const ulong CUTOFF5 =   84;
     603     2624828 :   const ulong CUTOFF6 =   40;
     604     2624828 :   const ulong CUTOFF7 =   23;
     605     2624828 :   const ulong CUTOFF8 =   15;
     606     2624828 :   const ulong CUTOFF9 =   11;
     607     2624828 :   const ulong CUTOFF10 =   9;
     608     2624828 :   const ulong CUTOFF11 =   7;
     609     2624828 :   const ulong CUTOFF12 =   6;
     610     2624828 :   const ulong CUTOFF13 =   5;
     611     2624828 :   const ulong CUTOFF14 =   4;
     612     2624828 :   const ulong CUTOFF15 =   4;
     613     2624828 :   const ulong CUTOFF16 =   3;
     614     2624828 :   const ulong CUTOFF17 =   3;
     615     2624828 :   const ulong CUTOFF18 =   3;
     616     2624828 :   const ulong CUTOFF19 =   3;
     617     2624828 :   const ulong CUTOFF20 =   3;
     618             : #endif
     619             : 
     620    18256218 :   if (p <= 2)
     621             :   {
     622      778233 :     if (p < 2) return p;
     623      775506 :     return k < BITS_IN_LONG? 1UL<<k: 0;
     624             :   }
     625    17477985 :   switch(k)
     626             :   {
     627             :     ulong p2, p3, p4, p5, p8;
     628     1731044 :     case 0:  return 1;
     629     7528237 :     case 1:  return p;
     630     2716996 :     case 2:  return usqru(p);
     631     1628329 :     case 3:  if (p > CUTOFF3) return 0; return p*p*p;
     632      829392 :     case 4:  if (p > CUTOFF4) return 0; p2=p*p; return p2*p2;
     633      464517 :     case 5:  if (p > CUTOFF5) return 0; p2=p*p; return p2*p2*p;
     634      551517 :     case 6:  if (p > CUTOFF6) return 0; p2=p*p; return p2*p2*p2;
     635      140761 :     case 7:  if (p > CUTOFF7) return 0; p2=p*p; return p2*p2*p2*p;
     636      120730 :     case 8:  if (p > CUTOFF8) return 0; p2=p*p; p4=p2*p2; return p4*p4;
     637      207911 :     case 9:  if (p > CUTOFF9) return 0; p2=p*p; p4=p2*p2; return p4*p4*p;
     638       68231 :     case 10: if (p > CUTOFF10)return 0; p2=p*p; p4=p2*p2; return p4*p4*p2;
     639       34215 :     case 11: if (p > CUTOFF11)return 0; p2=p*p; p4=p2*p2; return p4*p4*p2*p;
     640       51837 :     case 12: if (p > CUTOFF12)return 0; p2=p*p; p4=p2*p2; return p4*p4*p4;
     641       61926 :     case 13: if (p > CUTOFF13)return 0; p2=p*p; p4=p2*p2; return p4*p4*p4*p;
     642       66963 :     case 14: if (p > CUTOFF14)return 0; p2=p*p; p4=p2*p2; return p4*p4*p4*p2;
     643       86329 :     case 15: if (p > CUTOFF15)return 0;
     644       18914 :       p2=p*p; p3=p2*p; p5=p3*p2; return p5*p5*p5;
     645       60328 :     case 16: if (p > CUTOFF16)return 0;
     646       14808 :       p2=p*p; p4=p2*p2; p8=p4*p4; return p8*p8;
     647       34709 :     case 17: if (p > CUTOFF17)return 0;
     648        9721 :       p2=p*p; p4=p2*p2; p8=p4*p4; return p*p8*p8;
     649       24334 :     case 18: if (p > CUTOFF18)return 0;
     650        7254 :       p2=p*p; p4=p2*p2; p8=p4*p4; return p2*p8*p8;
     651      648272 :     case 19: if (p > CUTOFF19)return 0;
     652      634829 :       p2=p*p; p4=p2*p2; p8=p4*p4; return p*p2*p8*p8;
     653       12671 :     case 20: if (p > CUTOFF20)return 0;
     654        3905 :       p2=p*p; p4=p2*p2; p8=p4*p4; return p4*p8*p8;
     655             :   }
     656             : #ifdef LONG_IS_64BIT
     657      349687 :   switch(p)
     658             :   {
     659       35556 :     case 3: if (k > 40) return 0;
     660       21432 :       break;
     661         876 :     case 4: if (k > 31) return 0;
     662          90 :       return 1UL<<(2*k);
     663       15535 :     case 5: if (k > 27) return 0;
     664        3007 :       break;
     665         762 :     case 6: if (k > 24) return 0;
     666          42 :       break;
     667       10308 :     case 7: if (k > 22) return 0;
     668        1044 :       break;
     669      286650 :     default: return 0;
     670             :   }
     671             :   /* no overflow */
     672             :   {
     673       25525 :     ulong q = upowuu(p, k >> 1);
     674       25525 :     q *= q ;
     675       25525 :     return odd(k)? q*p: q;
     676             :   }
     677             : #else
     678       59049 :   return 0;
     679             : #endif
     680             : }
     681             : 
     682             : typedef struct {
     683             :   long prec, a;
     684             :   GEN (*sqr)(GEN);
     685             :   GEN (*mulug)(ulong,GEN);
     686             : } sr_muldata;
     687             : 
     688             : static GEN
     689      465821 : _rpowuu_sqr(void *data, GEN x)
     690             : {
     691      465821 :   sr_muldata *D = (sr_muldata *)data;
     692      465821 :   if (typ(x) == t_INT && lgefint(x) >= D->prec)
     693             :   { /* switch to t_REAL */
     694       18759 :     D->sqr   = &sqrr;
     695       18759 :     D->mulug = &mulur; x = itor(x, D->prec);
     696             :   }
     697      465821 :   return D->sqr(x);
     698             : }
     699             : 
     700             : static GEN
     701      178236 : _rpowuu_msqr(void *data, GEN x)
     702             : {
     703      178236 :   GEN x2 = _rpowuu_sqr(data, x);
     704      178236 :   sr_muldata *D = (sr_muldata *)data;
     705      178236 :   return D->mulug(D->a, x2);
     706             : }
     707             : 
     708             : /* return a^n as a t_REAL of precision prec. Assume a > 0, n > 0 */
     709             : GEN
     710      125519 : rpowuu(ulong a, ulong n, long prec)
     711             : {
     712             :   pari_sp av;
     713             :   GEN y, z;
     714             :   sr_muldata D;
     715             : 
     716      125519 :   if (a == 1) return real_1(prec);
     717      125519 :   if (a == 2) return real2n(n, prec);
     718      125519 :   if (n == 1) return utor(a, prec);
     719      124714 :   z = cgetr(prec);
     720      124714 :   av = avma;
     721      124714 :   D.sqr   = &sqri;
     722      124714 :   D.mulug = &mului;
     723      124714 :   D.prec = prec;
     724      124714 :   D.a = (long)a;
     725      124714 :   y = gen_powu_fold_i(utoipos(a), n, (void*)&D, &_rpowuu_sqr, &_rpowuu_msqr);
     726      124714 :   mpaff(y, z); avma = av; return z;
     727             : }
     728             : 
     729             : GEN
     730     5389104 : powrs(GEN x, long n)
     731             : {
     732     5389104 :   pari_sp av = avma;
     733             :   GEN y;
     734     5389104 :   if (!n) return powr0(x);
     735     5389104 :   y = gen_powu_i(x, (ulong)labs(n), NULL, &_sqrr, &_mulr);
     736     5389104 :   if (n < 0) y = invr(y);
     737     5389104 :   return gerepileuptoleaf(av,y);
     738             : }
     739             : GEN
     740      544932 : powru(GEN x, ulong n)
     741             : {
     742      544932 :   pari_sp av = avma;
     743             :   GEN y;
     744      544932 :   if (!n) return powr0(x);
     745      543735 :   y = gen_powu_i(x, n, NULL, &_sqrr, &_mulr);
     746      543735 :   return gerepileuptoleaf(av,y);
     747             : }
     748             : 
     749             : GEN
     750       24011 : powersr(GEN x, long n)
     751             : {
     752       24011 :   long prec = realprec(x);
     753       24011 :   return gen_powers(x, n, 1, &prec, &_sqrr, &_mulr, &_oner);
     754             : }
     755             : 
     756             : /* x^(s/2), assume x t_REAL */
     757             : GEN
     758           0 : powrshalf(GEN x, long s)
     759             : {
     760           0 :   if (s & 1) return sqrtr(powrs(x, s));
     761           0 :   return powrs(x, s>>1);
     762             : }
     763             : /* x^(s/2), assume x t_REAL */
     764             : GEN
     765       18915 : powruhalf(GEN x, ulong s)
     766             : {
     767       18915 :   if (s & 1) return sqrtr(powru(x, s));
     768        5439 :   return powru(x, s>>1);
     769             : }
     770             : /* x^(n/d), assume x t_REAL, return t_REAL */
     771             : GEN
     772         504 : powrfrac(GEN x, long n, long d)
     773             : {
     774             :   long z;
     775         504 :   if (!n) return powr0(x);
     776           0 :   z = cgcd(n, d); if (z > 1) { n /= z; d /= z; }
     777           0 :   if (d == 1) return powrs(x, n);
     778           0 :   x = powrs(x, n);
     779           0 :   if (d == 2) return sqrtr(x);
     780           0 :   return sqrtnr(x, d);
     781             : }
     782             : 
     783             : /* assume x != 0 */
     784             : static GEN
     785      559039 : pow_monome(GEN x, long n)
     786             : {
     787      559039 :   long i, d, dx = degpol(x);
     788             :   GEN A, b, y;
     789             : 
     790      559039 :   if (n < 0) { n = -n; y = cgetg(3, t_RFRAC); } else y = NULL;
     791             : 
     792      559039 :   if (HIGHWORD(dx) || HIGHWORD(n))
     793           8 :   {
     794             :     LOCAL_HIREMAINDER;
     795           9 :     d = (long)mulll((ulong)dx, (ulong)n);
     796           9 :     if (hiremainder || (d &~ LGBITS)) d = LGBITS; /* overflow */
     797           9 :     d += 2;
     798             :   }
     799             :   else
     800      559030 :     d = dx*n + 2;
     801      559039 :   if ((d + 1) & ~LGBITS) pari_err(e_OVERFLOW,"pow_monome [degree]");
     802      559032 :   A = cgetg(d+1, t_POL); A[1] = x[1];
     803      559032 :   for (i=2; i < d; i++) gel(A,i) = gen_0;
     804      559032 :   b = gpowgs(gel(x,dx+2), n); /* not memory clean if (n < 0) */
     805      559032 :   if (!y) y = A;
     806             :   else {
     807       20433 :     GEN c = denom(b);
     808       20433 :     gel(y,1) = c; if (c != gen_1) b = gmul(b,c);
     809       20433 :     gel(y,2) = A;
     810             :   }
     811      559032 :   gel(A,d) = b; return y;
     812             : }
     813             : 
     814             : /* x t_PADIC */
     815             : static GEN
     816       40971 : powps(GEN x, long n)
     817             : {
     818       40971 :   long e = n*valp(x), v;
     819       40971 :   GEN t, y, mod, p = gel(x,2);
     820             :   pari_sp av;
     821             : 
     822       40971 :   if (!signe(gel(x,4))) {
     823          84 :     if (n < 0) pari_err_INV("powps",x);
     824          77 :     return zeropadic(p, e);
     825             :   }
     826       40887 :   v = z_pval(n, p);
     827             : 
     828       40887 :   y = cgetg(5,t_PADIC);
     829       40887 :   mod = gel(x,3);
     830       40887 :   if (v == 0) mod = icopy(mod);
     831             :   else
     832             :   {
     833       39837 :     if (precp(x) == 1 && absequaliu(p, 2)) v++;
     834       39837 :     mod = mulii(mod, powiu(p,v));
     835       39837 :     mod = gerepileuptoint((pari_sp)y, mod);
     836             :   }
     837       40887 :   y[1] = evalprecp(precp(x) + v) | evalvalp(e);
     838       40887 :   gel(y,2) = icopy(p);
     839       40887 :   gel(y,3) = mod;
     840             : 
     841       40887 :   av = avma; t = gel(x,4);
     842       40887 :   if (n < 0) { t = Fp_inv(t, mod); n = -n; }
     843       40887 :   t = Fp_powu(t, n, mod);
     844       40887 :   gel(y,4) = gerepileuptoint(av, t);
     845       40887 :   return y;
     846             : }
     847             : /* x t_PADIC */
     848             : static GEN
     849         161 : powp(GEN x, GEN n)
     850             : {
     851             :   long v;
     852         161 :   GEN y, mod, p = gel(x,2);
     853             : 
     854         161 :   if (valp(x)) pari_err_OVERFLOW("valp()");
     855             : 
     856         161 :   if (!signe(gel(x,4))) {
     857          14 :     if (signe(n) < 0) pari_err_INV("powp",x);
     858           7 :     return zeropadic(p, 0);
     859             :   }
     860         147 :   v = Z_pval(n, p);
     861             : 
     862         147 :   y = cgetg(5,t_PADIC);
     863         147 :   mod = gel(x,3);
     864         147 :   if (v == 0) mod = icopy(mod);
     865             :   else
     866             :   {
     867          70 :     mod = mulii(mod, powiu(p,v));
     868          70 :     mod = gerepileuptoint((pari_sp)y, mod);
     869             :   }
     870         147 :   y[1] = evalprecp(precp(x) + v) | _evalvalp(0);
     871         147 :   gel(y,2) = icopy(p);
     872         147 :   gel(y,3) = mod;
     873         147 :   gel(y,4) = Fp_pow(gel(x,4), n, mod);
     874         147 :   return y;
     875             : }
     876             : static GEN
     877       40193 : pow_polmod(GEN x, GEN n)
     878             : {
     879       40193 :   GEN z = cgetg(3, t_POLMOD), a = gel(x,2), T = gel(x,1);
     880       40193 :   gel(z,1) = gcopy(T);
     881       40193 :   if (typ(a) != t_POL || varn(a) != varn(T) || lg(a) <= 3)
     882        2541 :     a = powgi(a, n);
     883             :   else {
     884       37652 :     pari_sp av = avma;
     885       37652 :     GEN p = NULL;
     886       37652 :     if (RgX_is_FpX(T, &p) && RgX_is_FpX(a, &p) && p)
     887             :     {
     888        7602 :       T = RgX_to_FpX(T, p); a = RgX_to_FpX(a, p);
     889        7602 :       if (lgefint(p) == 3)
     890             :       {
     891        7595 :         ulong pp = p[2];
     892        7595 :         a = Flxq_pow(ZX_to_Flx(a, pp), n, ZX_to_Flx(T, pp), pp);
     893        7595 :         a = Flx_to_ZX(a);
     894             :       }
     895             :       else
     896           7 :         a = FpXQ_pow(a, n, T, p);
     897        7602 :       a = FpX_to_mod(a, p);
     898        7602 :       a = gerepileupto(av, a);
     899             :     }
     900             :     else
     901             :     {
     902       30050 :       avma = av;
     903       30050 :       a = RgXQ_pow(a, n, gel(z,1));
     904             :     }
     905             :   }
     906       40193 :   gel(z,2) = a; return z;
     907             : }
     908             : 
     909             : GEN
     910    45636706 : gpowgs(GEN x, long n)
     911             : {
     912             :   long m;
     913             :   pari_sp av;
     914             :   GEN y;
     915             : 
     916    45636706 :   if (n == 0) return gpowg0(x);
     917    42863689 :   if (n == 1)
     918     4190016 :     switch (typ(x)) {
     919      673395 :       case t_QFI: return redimag(x);
     920          14 :       case t_QFR: return redreal(x);
     921     3516607 :       default: return gcopy(x);
     922             :     }
     923    38673673 :   if (n ==-1) return ginv(x);
     924    35276801 :   switch(typ(x))
     925             :   {
     926    15554205 :     case t_INT: return powis(x,n);
     927     5388593 :     case t_REAL: return powrs(x,n);
     928             :     case t_INTMOD:
     929       21000 :       y = cgetg(3,t_INTMOD); gel(y,1) = icopy(gel(x,1));
     930       21000 :       gel(y,2) = Fp_pows(gel(x,2), n, gel(x,1));
     931       21000 :       return y;
     932             :     case t_FRAC:
     933             :     {
     934      224135 :       GEN a = gel(x,1), b = gel(x,2);
     935      224135 :       long s = (signe(a) < 0 && odd(n))? -1: 1;
     936      224135 :       if (n < 0) {
     937          14 :         n = -n;
     938          14 :         if (is_pm1(a)) return powiu_sign(b, n, s); /* +-1/x[2] inverts to t_INT */
     939           7 :         swap(a, b);
     940             :       }
     941      224128 :       y = cgetg(3, t_FRAC);
     942      224128 :       gel(y,1) = powiu_sign(a, n, s);
     943      224128 :       gel(y,2) = powiu_sign(b, n, 1);
     944      224128 :       return y;
     945             :     }
     946       40971 :     case t_PADIC: return powps(x, n);
     947             :     case t_RFRAC:
     948             :     {
     949      195685 :       av = avma; y = cgetg(3, t_RFRAC); m = labs(n);
     950      195685 :       gel(y,1) = gpowgs(gel(x,1),m);
     951      195685 :       gel(y,2) = gpowgs(gel(x,2),m);
     952      195685 :       if (n < 0) y = ginv(y);
     953      195685 :       return gerepileupto(av,y);
     954             :     }
     955             :     case t_POLMOD: {
     956       40186 :       long N[] = {evaltyp(t_INT) | _evallg(3),0,0};
     957       40186 :       affsi(n,N); return pow_polmod(x, N);
     958             :     }
     959             :     case t_POL:
     960      622060 :       if (RgX_is_monomial(x)) return pow_monome(x, n);
     961             :     default: {
     962    13252987 :       pari_sp av = avma;
     963    13252987 :       y = gen_powu_i(x, (ulong)labs(n), NULL, &_sqr, &_mul);
     964    13252986 :       if (n < 0) y = ginv(y);
     965    13252986 :       return gerepileupto(av,y);
     966             :     }
     967             :   }
     968             : }
     969             : 
     970             : /* n a t_INT */
     971             : GEN
     972    33199312 : powgi(GEN x, GEN n)
     973             : {
     974             :   GEN y;
     975             : 
     976    33199312 :   if (!is_bigint(n)) return gpowgs(x, itos(n));
     977             :   /* probable overflow for non-modular types (typical exception: (X^0)^N) */
     978         340 :   switch(typ(x))
     979             :   {
     980             :     case t_INTMOD:
     981          49 :       y = cgetg(3,t_INTMOD); gel(y,1) = icopy(gel(x,1));
     982          49 :       gel(y,2) = Fp_pow(gel(x,2), n, gel(x,1));
     983          49 :       return y;
     984          58 :     case t_FFELT: return FF_pow(x,n);
     985         161 :     case t_PADIC: return powp(x, n);
     986             : 
     987             :     case t_INT:
     988          35 :       if (is_pm1(x)) return (signe(x) < 0 && mpodd(n))? gen_m1: gen_1;
     989          14 :       if (signe(x)) pari_err_OVERFLOW("lg()");
     990           7 :       if (signe(n) < 0) pari_err_INV("powgi",gen_0);
     991           7 :       return gen_0;
     992             :     case t_FRAC:
     993           7 :       pari_err_OVERFLOW("lg()");
     994             : 
     995          14 :     case t_QFR: return qfrpow(x, n);
     996           7 :     case t_POLMOD: return pow_polmod(x, n);
     997             :     default: {
     998           9 :       pari_sp av = avma;
     999           9 :       y = gen_pow(x, n, NULL, &_sqr, &_mul);
    1000           9 :       if (signe(n) < 0) y = ginv(y);
    1001           9 :       return gerepileupto(av,y);
    1002             :     }
    1003             :   }
    1004             : }
    1005             : 
    1006             : /* Assume x = 1 + O(t), n a scalar. Return x^n */
    1007             : static GEN
    1008         196 : ser_pow_1(GEN x, GEN n)
    1009             : {
    1010             :   long lx, mi, i, j, d;
    1011         196 :   GEN y = cgetg_copy(x, &lx), X = x+2, Y = y + 2;
    1012         196 :   y[1] = evalsigne(1) | _evalvalp(0) | evalvarn(varn(x));
    1013         196 :   d = mi = lx-3; while (mi>=1 && isrationalzero(gel(X,mi))) mi--;
    1014         196 :   gel(Y,0) = gen_1;
    1015        1974 :   for (i=1; i<=d; i++)
    1016             :   {
    1017        1778 :     pari_sp av = avma;
    1018        1778 :     GEN s = gen_0;
    1019       15246 :     for (j=1; j<=minss(i,mi); j++)
    1020             :     {
    1021       13468 :       GEN t = gsubgs(gmulgs(n,j),i-j);
    1022       13468 :       s = gadd(s, gmul(gmul(t, gel(X,j)), gel(Y,i-j)));
    1023             :     }
    1024        1778 :     gel(Y,i) = gerepileupto(av, gdivgs(s,i));
    1025             :   }
    1026         196 :   return y;
    1027             : }
    1028             : 
    1029             : /* we suppose n != 0, valp(x) = 0 and leading-term(x) != 0. Not stack clean */
    1030             : static GEN
    1031         196 : ser_pow(GEN x, GEN n, long prec)
    1032             : {
    1033             :   GEN y, c, lead;
    1034         196 :   if (varncmp(gvar(n), varn(x)) <= 0) return gexp(gmul(n, glog(x,prec)), prec);
    1035         196 :   lead = gel(x,2);
    1036         196 :   if (gequal1(lead)) return ser_pow_1(x, n);
    1037          91 :   x = ser_normalize(x);
    1038          91 :   if (typ(n) == t_FRAC && !isinexact(lead) && ispower(lead, gel(n,2), &c))
    1039          84 :     c = powgi(c, gel(n,1));
    1040             :   else
    1041           7 :     c = gpow(lead,n, prec);
    1042          91 :   y = gmul(c, ser_pow_1(x, n));
    1043             :   /* gpow(t_POLMOD,n) can be a t_COL [conjvec] */
    1044          91 :   if (typ(y) != t_SER) pari_err_TYPE("gpow", y);
    1045          91 :   return y;
    1046             : }
    1047             : 
    1048             : static long
    1049         210 : val_from_i(GEN E)
    1050             : {
    1051         210 :   if (is_bigint(E)) pari_err_OVERFLOW("sqrtn [valuation]");
    1052         203 :   return itos(E);
    1053             : }
    1054             : 
    1055             : /* return x^q, assume typ(x) = t_SER, typ(q) = t_INT/t_FRAC and q != 0 */
    1056             : static GEN
    1057         217 : ser_powfrac(GEN x, GEN q, long prec)
    1058             : {
    1059         217 :   GEN y, E = gmulsg(valp(x), q);
    1060             :   long e;
    1061             : 
    1062         217 :   if (!signe(x))
    1063             :   {
    1064          21 :     if (gsigne(q) < 0) pari_err_INV("gpow", x);
    1065          21 :     return zeroser(varn(x), val_from_i(gfloor(E)));
    1066             :   }
    1067         196 :   if (typ(E) != t_INT)
    1068           7 :     pari_err_DOMAIN("sqrtn", "valuation", "!=", mkintmod(gen_0, gel(q,2)), x);
    1069         189 :   e = val_from_i(E);
    1070         189 :   y = leafcopy(x); setvalp(y, 0);
    1071         189 :   y = ser_pow(y, q, prec);
    1072         189 :   setvalp(y, e); return y;
    1073             : }
    1074             : 
    1075             : static GEN
    1076         126 : gpow0(GEN x, GEN n, long prec)
    1077             : {
    1078         126 :   pari_sp av = avma;
    1079             :   long i, lx;
    1080             :   GEN y;
    1081         126 :   switch(typ(n))
    1082             :   {
    1083             :     case t_INT: case t_REAL: case t_FRAC: case t_COMPLEX: case t_QUAD:
    1084          84 :       break;
    1085             :     case t_VEC: case t_COL: case t_MAT:
    1086          35 :       y = cgetg_copy(n, &lx);
    1087          35 :       for (i=1; i<lx; i++) gel(y,i) = gpow0(x,gel(n,i),prec);
    1088          35 :       return y;
    1089           7 :     default: pari_err_TYPE("gpow(0,n)", n);
    1090             :   }
    1091          84 :   n = real_i(n);
    1092          84 :   if (gsigne(n) <= 0) pari_err_DOMAIN("gpow(0,n)", "n", "<=", gen_0, n);
    1093          77 :   if (!precision(x)) return gcopy(x);
    1094             : 
    1095          14 :   x = ground(gmulsg(gexpo(x),n));
    1096          14 :   if (is_bigint(x) || uel(x,2) >= HIGHEXPOBIT)
    1097           7 :     pari_err_OVERFLOW("gpow");
    1098           7 :   avma = av; return real_0_bit(itos(x));
    1099             : }
    1100             : 
    1101             : GEN
    1102    12826916 : gpow(GEN x, GEN n, long prec)
    1103             : {
    1104    12826916 :   long i, lx, tx, tn = typ(n);
    1105             :   pari_sp av;
    1106             :   GEN y;
    1107             : 
    1108    12826916 :   if (tn == t_INT) return powgi(x,n);
    1109     2159990 :   tx = typ(x);
    1110     2159990 :   if (is_matvec_t(tx))
    1111             :   {
    1112          49 :     y = cgetg_copy(x, &lx);
    1113          49 :     for (i=1; i<lx; i++) gel(y,i) = gpow(gel(x,i),n,prec);
    1114          49 :     return y;
    1115             :   }
    1116     2159941 :   av = avma;
    1117     2159941 :   switch (tx)
    1118             :   {
    1119          14 :     case t_POL: case t_RFRAC: x = toser_i(x); /* fall through */
    1120             :     case t_SER:
    1121          77 :       if (tn == t_FRAC) return gerepileupto(av, ser_powfrac(x, n, prec));
    1122          35 :       if (valp(x))
    1123          21 :         pari_err_DOMAIN("gpow [irrational exponent]",
    1124             :                         "valuation", "!=", gen_0, x);
    1125          14 :       if (lg(x) == 2) return gerepilecopy(av, x); /* O(1) */
    1126           7 :       return gerepileupto(av, ser_pow(x, n, prec));
    1127             :   }
    1128     2159864 :   if (gequal0(x)) return gpow0(x, n, prec);
    1129     2159808 :   if (tn == t_FRAC)
    1130             :   {
    1131     1942319 :     GEN z, d = gel(n,2), a = gel(n,1);
    1132             :     long D;
    1133     1942319 :     switch (tx)
    1134             :     {
    1135             :     case t_INTMOD:
    1136             :       {
    1137          21 :         GEN p = gel(x,1);
    1138          21 :         if (!BPSW_psp(p)) pari_err_PRIME("gpow",p);
    1139          14 :         y = cgetg(3,t_INTMOD); gel(y,1) = icopy(p);
    1140          14 :         av = avma;
    1141          14 :         z = Fp_sqrtn(gel(x,2), d, p, NULL);
    1142          14 :         if (!z) pari_err_SQRTN("gpow",x);
    1143           7 :         gel(y,2) = gerepileuptoint(av, Fp_pow(z, a, p));
    1144           7 :         return y;
    1145             :       }
    1146             : 
    1147             :     case t_PADIC:
    1148          14 :       z = Qp_sqrtn(x, d, NULL); if (!z) pari_err_SQRTN("gpow",x);
    1149           7 :       return gerepileupto(av, powgi(z, a));
    1150             : 
    1151             :     case t_FFELT:
    1152          21 :       return gerepileupto(av,FF_pow(FF_sqrtn(x,d,NULL),a));
    1153             :     default:
    1154     1942263 :       D = itos_or_0(d);
    1155     1942263 :       if (!D) break;
    1156     1942263 :       if (D == 2)
    1157             :       {
    1158     1685997 :         GEN y = gsqrt(x,prec), t = shifti(subiu(a,1), -1);
    1159     1685997 :         if (signe(t)) y = gmul(y, powgi(x,t));
    1160     1685997 :         return gerepileupto(av, y);
    1161             :       }
    1162      256266 :       if (is_real_t(tx) && gsigne(x) > 0)
    1163             :       {
    1164      255395 :         if (tx != t_REAL) x = gtofp(x, prec);
    1165      255395 :         z = sqrtnr(x, D);
    1166      255395 :         if (!equali1(a)) z = powgi(z, a);
    1167      255395 :         return gerepileuptoleaf(av, z);
    1168             :       }
    1169             :     }
    1170             :   }
    1171      218360 :   i = (long) precision(n); if (i) prec=i;
    1172      218360 :   y = gmul(n, glog(x,prec));
    1173      218304 :   return gerepileupto(av, gexp(y,prec));
    1174             : }
    1175             : 
    1176             : GEN
    1177       21444 : gpowers0(GEN x, long n, GEN x0)
    1178             : {
    1179             :   long i, l;
    1180             :   GEN V;
    1181       21444 :   if (!x0) return gpowers(x,n);
    1182       21416 :   if (n < 0) return cgetg(1,t_VEC);
    1183       21416 :   l = n+2; V = cgetg(l, t_VEC); gel(V,1) = gcopy(x0);
    1184       21416 :   for (i = 2; i < l; i++) gel(V,i) = gmul(gel(V,i-1),x);
    1185       21416 :   return V;
    1186             : }
    1187             : 
    1188             : GEN
    1189       45853 : gpowers(GEN x, long n)
    1190             : {
    1191       45853 :   if (n < 0) return cgetg(1,t_VEC);
    1192       45846 :   return gen_powers(x, n, 1, (void*)x, &_sqr, &_mul, &_one);
    1193             : }
    1194             : 
    1195             : /* return [q^1,q^4,...,q^{n^2}] */
    1196             : GEN
    1197       19085 : gsqrpowers(GEN q, long n)
    1198             : {
    1199       19085 :   pari_sp av = avma;
    1200       19085 :   GEN L = gpowers0(gsqr(q), n, q); /* L[i] = q^(2i - 1), i <= n+1 */
    1201       19085 :   GEN v = cgetg(n+1, t_VEC);
    1202             :   long i;
    1203       19085 :   gel(v, 1) = gcopy(q);
    1204       19085 :   for (i = 2; i <= n ; ++i) gel(v, i) = q = gmul(q, gel(L,i)); /* q^(i^2) */
    1205       19085 :   return gerepileupto(av, v);
    1206             : }
    1207             : 
    1208             : /* 4 | N. returns a vector RU which contains exp(2*i*k*Pi/N), k=0..N-1 */
    1209             : static GEN
    1210       56479 : grootsof1_4(long N, long prec)
    1211             : {
    1212       56479 :   GEN z, RU = cgetg(N+1,t_VEC), *v  = ((GEN*)RU) + 1;
    1213       56479 :   long i, N2 = (N>>1), N4 = (N>>2), N8 = (N>>3);
    1214             :   /* z^N2 = -1, z^N4 = I; if z^k = a+I*b, then z^(N4-k) = I*conj(z) = b+a*I */
    1215             : 
    1216       56479 :   v[0] = gen_1; v[1] = z = rootsof1u_cx(N, prec);
    1217       56479 :   if (odd(N4)) N8++;
    1218       78057 :   for (i=1; i<N8; i++)
    1219             :   {
    1220       21578 :     GEN t = v[i];
    1221       21578 :     v[i+1] = gmul(z, t);
    1222       21578 :     v[N4-i] = mkcomplex(gel(t,2), gel(t,1));
    1223             :   }
    1224       56479 :   for (i=0; i<N4; i++) v[i+N4] = mulcxI(v[i]);
    1225       56479 :   for (i=0; i<N2; i++) v[i+N2] = gneg(v[i]);
    1226       56479 :   return RU;
    1227             : }
    1228             : 
    1229             : /* as above, N arbitrary */
    1230             : GEN
    1231       62656 : grootsof1(long N, long prec)
    1232             : {
    1233             :   GEN z, RU, *v;
    1234             :   long i, k;
    1235             : 
    1236       62656 :   if ((N & 3) == 0) return grootsof1_4(N, prec);
    1237        6177 :   if (N == 1) return mkvec(gen_1);
    1238        6170 :   k = (N+3)>>1;
    1239        6170 :   RU = cgetg(N+1,t_VEC);
    1240        6170 :   v  = ((GEN*)RU) + 1;
    1241        6170 :   v[0] = gen_1; v[1] = z = rootsof1u_cx(N, prec);
    1242        6170 :   for (i=2; i<k; i++) v[i] = gmul(z, v[i-1]);
    1243        6170 :   for (   ; i<N; i++) v[i] = gconj(v[N-i]);
    1244        6170 :   return RU;
    1245             : }
    1246             : 
    1247             : /********************************************************************/
    1248             : /**                                                                **/
    1249             : /**                        RACINE CARREE                           **/
    1250             : /**                                                                **/
    1251             : /********************************************************************/
    1252             : /* assume x unit, e = precp(x) */
    1253             : GEN
    1254        1477 : Z2_sqrt(GEN x, long e)
    1255             : {
    1256        1477 :   ulong r = signe(x)>=0?mod16(x):16-mod16(x);
    1257             :   GEN z;
    1258             :   long ez;
    1259             :   pari_sp av;
    1260             : 
    1261        1477 :   switch(e)
    1262             :   {
    1263           7 :     case 1: return gen_1;
    1264         119 :     case 2: return (r & 3UL) == 1? gen_1: NULL;
    1265          14 :     case 3: return (r & 7UL) == 1? gen_1: NULL;
    1266          70 :     case 4: if (r == 1) return gen_1;
    1267          35 :             else return (r == 9)? utoipos(3): NULL;
    1268        1267 :     default: if ((r&7UL) != 1) return NULL;
    1269             :   }
    1270        1267 :   av = avma; z = (r==1)? gen_1: utoipos(3);
    1271        1267 :   ez = 3; /* number of correct bits in z (compared to sqrt(x)) */
    1272             :   for(;;)
    1273             :   {
    1274             :     GEN mod;
    1275        4634 :     ez = (ez<<1) - 1;
    1276        4634 :     if (ez > e) ez = e;
    1277        4634 :     mod = int2n(ez);
    1278        4634 :     z = addii(z, remi2n(mulii(x, Fp_inv(z,mod)), ez));
    1279        4634 :     z = shifti(z, -1); /* (z + x/z) / 2 */
    1280        4634 :     if (e == ez) return gerepileuptoint(av, z);
    1281        3367 :     if (ez < e) ez--;
    1282        3367 :     if (gc_needed(av,2))
    1283             :     {
    1284           0 :       if (DEBUGMEM > 1) pari_warn(warnmem,"Qp_sqrt");
    1285           0 :       z = gerepileuptoint(av,z);
    1286             :     }
    1287        3367 :   }
    1288             : }
    1289             : 
    1290             : /* x unit defined modulo p^e, e > 0 */
    1291             : GEN
    1292        1736 : Qp_sqrt(GEN x)
    1293             : {
    1294        1736 :   long pp, e = valp(x);
    1295        1736 :   GEN z,y,mod, p = gel(x,2);
    1296             : 
    1297        1736 :   if (gequal0(x)) return zeropadic(p, (e+1) >> 1);
    1298        1736 :   if (e & 1) return NULL;
    1299             : 
    1300        1722 :   y = cgetg(5,t_PADIC);
    1301        1722 :   pp = precp(x);
    1302        1722 :   mod = gel(x,3);
    1303        1722 :   z   = gel(x,4); /* lift to t_INT */
    1304        1722 :   e >>= 1;
    1305        1722 :   z = Zp_sqrt(z, p, pp);
    1306        1722 :   if (!z) return NULL;
    1307        1673 :   if (absequaliu(p,2))
    1308             :   {
    1309         805 :     pp  = (pp <= 3) ? 1 : pp-1;
    1310         805 :     mod = int2n(pp);
    1311             :   }
    1312         868 :   else mod = icopy(mod);
    1313        1673 :   y[1] = evalprecp(pp) | evalvalp(e);
    1314        1673 :   gel(y,2) = icopy(p);
    1315        1673 :   gel(y,3) = mod;
    1316        1673 :   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        1421 : sqrt_ser(GEN b, long prec)
    1358             : {
    1359        1421 :   long e = valp(b), vx = varn(b), lx, lold, j;
    1360             :   ulong mask;
    1361             :   GEN a, x, lta, ltx;
    1362             : 
    1363        1421 :   if (!signe(b)) return zeroser(vx, e>>1);
    1364        1421 :   a = leafcopy(b);
    1365        1421 :   x = cgetg_copy(b, &lx);
    1366        1421 :   if (e & 1)
    1367           7 :     pari_err_DOMAIN("sqrtn", "valuation", "!=", mkintmod(gen_0, gen_2), x);
    1368        1414 :   a[1] = x[1] = evalsigne(1) | evalvarn(0) | _evalvalp(0);
    1369        1414 :   lta = gel(a,2);
    1370        1414 :   if (gequal1(lta)) ltx = lta;
    1371          91 :   else if (!issquareall(lta,&ltx)) ltx = gsqrt(lta,prec);
    1372        1414 :   gel(x,2) = ltx;
    1373        1414 :   for (j = 3; j < lx; j++) gel(x,j) = gen_0;
    1374        1414 :   setlg(x,3);
    1375        1414 :   mask = quadratic_prec_mask(lx - 2);
    1376        1414 :   lold = 1;
    1377        7472 :   while (mask > 1)
    1378             :   {
    1379        4644 :     GEN y, x2 = gmul2n(x,1);
    1380        4644 :     long l = lold << 1;
    1381             : 
    1382        4644 :     if (mask & 1) l--;
    1383        4644 :     mask >>= 1;
    1384        4644 :     setlg(a, l + 2);
    1385        4644 :     setlg(x, l + 2);
    1386        4644 :     y = sqr_ser_part(x, lold, l-1) - lold;
    1387        4644 :     for (j = lold+2; j < l+2; j++) gel(y,j) = gsub(gel(y,j), gel(a,j));
    1388        4644 :     y += lold; setvalp(y, lold);
    1389        4644 :     y = normalize(y);
    1390        4644 :     y = gsub(x, gdiv(y, x2)); /* = gmul2n(gsub(x, gdiv(a,x)), -1); */
    1391        4644 :     for (j = lold+2; j < l+2; j++) gel(x,j) = gel(y,j);
    1392        4644 :     lold = l;
    1393             :   }
    1394        1414 :   x[1] = evalsigne(1) | evalvarn(vx) | _evalvalp(e >> 1);
    1395        1414 :   return x;
    1396             : }
    1397             : 
    1398             : GEN
    1399     6879153 : gsqrt(GEN x, long prec)
    1400             : {
    1401             :   pari_sp av;
    1402             :   GEN y;
    1403             : 
    1404     6879153 :   switch(typ(x))
    1405             :   {
    1406             :     case t_INT:
    1407     1972330 :       if (!signe(x)) return real_0(prec); /* no loss of accuracy */
    1408     1972323 :       x = itor(x,prec); /* fall through */
    1409     6481773 :     case t_REAL: return sqrtr(x);
    1410             : 
    1411             :     case t_INTMOD:
    1412             :     {
    1413          35 :       GEN p = gel(x,1), a;
    1414          35 :       y = cgetg(3,t_INTMOD); gel(y,1) = icopy(p);
    1415          35 :       a = Fp_sqrt(gel(x,2),p);
    1416          21 :       if (!a)
    1417             :       {
    1418           7 :         if (!BPSW_psp(p)) pari_err_PRIME("sqrt [modulus]",p);
    1419           7 :         pari_err_SQRTN("gsqrt",x);
    1420             :       }
    1421          14 :       gel(y,2) = a; return y;
    1422             :     }
    1423             : 
    1424             :     case t_COMPLEX:
    1425             :     { /* (u+iv)^2 = a+ib <=> u^2+v^2 = sqrt(a^2+b^2), u^2-v^2=a, 2uv=b */
    1426      380683 :       GEN a = gel(x,1), b = gel(x,2), r, u, v;
    1427      380683 :       if (isrationalzero(b)) return gsqrt(a, prec);
    1428      380683 :       y = cgetg(3,t_COMPLEX); av = avma;
    1429             : 
    1430      380683 :       r = cxnorm(x);
    1431      380683 :       if (typ(r) == t_INTMOD) pari_err_IMPL("sqrt(complex of t_INTMODs)");
    1432      380683 :       r = gsqrt(r, prec); /* t_REAL, |a+Ib| */
    1433      380683 :       if (!signe(r))
    1434         117 :         u = v = gerepileuptoleaf(av, sqrtr(r));
    1435      380566 :       else if (gsigne(a) < 0)
    1436             :       {
    1437             :         /* v > 0 since r > 0, a < 0, rounding errors can't make the sum of two
    1438             :          * positive numbers = 0 */
    1439       39943 :         v = sqrtr( gmul2n(gsub(r,a), -1) );
    1440       39943 :         if (gsigne(b) < 0) togglesign(v);
    1441       39943 :         v = gerepileuptoleaf(av, v); av = avma;
    1442             :         /* v = 0 is impossible */
    1443       39943 :         u = gerepileuptoleaf(av, gdiv(b, shiftr(v,1)));
    1444             :       } else {
    1445      340623 :         u = sqrtr( gmul2n(gadd(r,a), -1) );
    1446      340623 :         u = gerepileuptoleaf(av, u); av = avma;
    1447      340623 :         if (!signe(u)) /* possible if a = 0.0, e.g. sqrt(0.e-10+1e-10*I) */
    1448           7 :           v = u;
    1449             :         else
    1450      340616 :           v = gerepileuptoleaf(av, gdiv(b, shiftr(u,1)));
    1451             :       }
    1452      380683 :       gel(y,1) = u;
    1453      380683 :       gel(y,2) = v; return y;
    1454             :     }
    1455             : 
    1456             :     case t_PADIC:
    1457          49 :       y = Qp_sqrt(x);
    1458          49 :       if (!y) pari_err_SQRTN("Qp_sqrt",x);
    1459          28 :       return y;
    1460             : 
    1461          70 :     case t_FFELT: return FF_sqrt(x);
    1462             : 
    1463             :     default:
    1464       16536 :       av = avma; if (!(y = toser_i(x))) break;
    1465        1421 :       return gerepilecopy(av, sqrt_ser(y, prec));
    1466             :   }
    1467       15115 :   return trans_eval("sqrt",gsqrt,x,prec);
    1468             : }
    1469             : /********************************************************************/
    1470             : /**                                                                **/
    1471             : /**                          N-th ROOT                             **/
    1472             : /**                                                                **/
    1473             : /********************************************************************/
    1474             : static void
    1475           7 : bug_logp(GEN p)
    1476             : {
    1477           7 :   if (!BPSW_psp(p)) pari_err_PRIME("p-adic log",p);
    1478           0 :   pari_err_BUG("log_p");
    1479           0 : }
    1480             : /* Let x = 1 mod p and y := (x-1)/(x+1) = 0 (p). Then
    1481             :  * log(x) = log(1+y) - log(1-y) = 2 \sum_{k odd} y^k / k.
    1482             :  * palogaux(x) returns the last sum (not multiplied by 2) */
    1483             : static GEN
    1484       23345 : palogaux(GEN x)
    1485             : {
    1486             :   long i, k, e, pp, t;
    1487       23345 :   GEN y,s,y2, p = gel(x,2);
    1488       23345 :   int is2 = absequaliu(p,2);
    1489             : 
    1490       23345 :   y = subiu(gel(x,4), 1);
    1491       23345 :   if (!signe(y))
    1492             :   {
    1493         812 :     long v = valp(x)+precp(x);
    1494         812 :     if (is2) v--;
    1495         812 :     return zeropadic(p, v);
    1496             :   }
    1497             :   /* optimize t: log(x) = log(x^(p^t)) / p^t */
    1498       22533 :   e = Z_pval(y, p); /* valp(y) = e >= 1; precp(y) = precp(x)-e */
    1499       22533 :   if (!e) bug_logp(p);
    1500       22526 :   if (is2)
    1501        2492 :     t = sqrt( (double)(precp(x)-e) / e ); /* instead of (2*e) */
    1502             :   else
    1503       20034 :     t = sqrt( (double)(precp(x)-e) / (e * (expi(p) + hammingweight(p))) );
    1504       22526 :   for (i = 0; i < t; i++) x = gpow(x, p, 0);
    1505             : 
    1506       22526 :   y = gdiv(gaddgs(x,-1), gaddgs(x,1));
    1507       22526 :   e = valp(y); /* > 0 */
    1508       22526 :   if (e <= 0) bug_logp(p);
    1509       22526 :   pp = precp(y) + e;
    1510       22526 :   if (is2) pp--;
    1511             :   else
    1512             :   {
    1513             :     GEN p1;
    1514       20034 :     for (p1=utoipos(e); abscmpui(pp,p1) > 0; pp++) p1 = mulii(p1, p);
    1515       20034 :     pp -= 2;
    1516             :   }
    1517       22526 :   k = pp/e; if (!odd(k)) k--;
    1518       22526 :   if (DEBUGLEVEL>5)
    1519           0 :     err_printf("logp: [pp,k,e,t] = [%ld,%ld,%ld,%ld]\n",pp,k,e,t);
    1520       22526 :   if (k > 1)
    1521             :   {
    1522       21763 :     y2 = gsqr(y); s = gdivgs(gen_1,k);
    1523      119930 :     while (k > 2)
    1524             :     {
    1525       76404 :       k -= 2;
    1526       76404 :       s = gadd(gmul(y2,s), gdivgs(gen_1,k));
    1527             :     }
    1528       21763 :     y = gmul(s,y);
    1529             :   }
    1530       22526 :   if (t) setvalp(y, valp(y) - t);
    1531       22526 :   return y;
    1532             : }
    1533             : 
    1534             : GEN
    1535       23352 : Qp_log(GEN x)
    1536             : {
    1537       23352 :   pari_sp av = avma;
    1538       23352 :   GEN y, p = gel(x,2), a = gel(x,4);
    1539             : 
    1540       23352 :   if (!signe(a)) pari_err_DOMAIN("Qp_log", "argument", "=", gen_0, x);
    1541       23345 :   y = leafcopy(x); setvalp(y,0);
    1542       23345 :   if (absequaliu(p,2))
    1543        2954 :     y = palogaux(gsqr(y));
    1544       20391 :   else if (gequal1(modii(a, p)))
    1545       11452 :     y = gmul2n(palogaux(y), 1);
    1546             :   else
    1547             :   { /* compute log(x^(p-1)) / (p-1) */
    1548        8939 :     GEN mod = gel(y,3), p1 = subiu(p,1);
    1549        8939 :     gel(y,4) = Fp_pow(a, p1, mod);
    1550        8939 :     p1 = diviiexact(subsi(1,mod), p1); /* 1/(p-1) */
    1551        8939 :     y = gmul(palogaux(y), shifti(p1,1));
    1552             :   }
    1553       23338 :   return gerepileupto(av,y);
    1554             : }
    1555             : 
    1556             : static GEN Qp_exp_safe(GEN x);
    1557             : 
    1558             : /*compute the p^e th root of x p-adic, assume x != 0 */
    1559             : static GEN
    1560         854 : Qp_sqrtn_ram(GEN x, long e)
    1561             : {
    1562         854 :   pari_sp ltop=avma;
    1563         854 :   GEN a, p = gel(x,2), n = powiu(p,e);
    1564         854 :   long v = valp(x), va;
    1565         854 :   if (v)
    1566             :   {
    1567             :     long z;
    1568         161 :     v = sdivsi_rem(v, n, &z);
    1569         161 :     if (z) return NULL;
    1570          91 :     x = leafcopy(x);
    1571          91 :     setvalp(x,0);
    1572             :   }
    1573             :   /*If p = 2, -1 is a root of 1 in U1: need extra check*/
    1574         784 :   if (absequaliu(p, 2) && mod8(gel(x,4)) != 1) return NULL;
    1575         749 :   a = Qp_log(x);
    1576         749 :   va = valp(a) - e;
    1577         749 :   if (va <= 0)
    1578             :   {
    1579         287 :     if (signe(gel(a,4))) return NULL;
    1580             :     /* all accuracy lost */
    1581         119 :     a = cvtop(remii(gel(x,4),p), p, 1);
    1582             :   }
    1583             :   else
    1584             :   {
    1585         462 :     setvalp(a, va); /* divide by p^e */
    1586         462 :     a = Qp_exp_safe(a);
    1587         462 :     if (!a) return NULL;
    1588             :     /* n=p^e and a^n=z*x where z is a (p-1)th-root of 1.
    1589             :      * Since z^n=z, we have (a/z)^n = x. */
    1590         462 :     a = gdiv(x, powgi(a,subiu(n,1))); /* = a/z = x/a^(n-1)*/
    1591         462 :     if (v) setvalp(a,v);
    1592             :   }
    1593         581 :   return gerepileupto(ltop,a);
    1594             : }
    1595             : 
    1596             : /*compute the nth root of x p-adic p prime with n*/
    1597             : static GEN
    1598         616 : Qp_sqrtn_unram(GEN x, GEN n, GEN *zetan)
    1599             : {
    1600             :   pari_sp av;
    1601         616 :   GEN Z, a, r, p = gel(x,2);
    1602         616 :   long v = valp(x);
    1603         616 :   if (v)
    1604             :   {
    1605             :     long z;
    1606          84 :     v = sdivsi_rem(v,n,&z);
    1607          84 :     if (z) return NULL;
    1608             :   }
    1609         609 :   r = cgetp(x); setvalp(r,v);
    1610         609 :   Z = NULL; /* -Wall */
    1611         609 :   if (zetan) Z = cgetp(x);
    1612         609 :   av = avma; a = Fp_sqrtn(gel(x,4), n, p, zetan);
    1613         609 :   if (!a) return NULL;
    1614         595 :   affii(Zp_sqrtnlift(gel(x,4), n, a, p, precp(x)), gel(r,4));
    1615         595 :   if (zetan)
    1616             :   {
    1617          14 :     affii(Zp_sqrtnlift(gen_1, n, *zetan, p, precp(x)), gel(Z,4));
    1618          14 :     *zetan = Z;
    1619             :   }
    1620         595 :   avma = av; return r;
    1621             : }
    1622             : 
    1623             : GEN
    1624        1183 : Qp_sqrtn(GEN x, GEN n, GEN *zetan)
    1625             : {
    1626             :   pari_sp av, tetpil;
    1627             :   GEN q, p;
    1628             :   long e;
    1629        1183 :   if (absequaliu(n, 2))
    1630             :   {
    1631          70 :     if (zetan) *zetan = gen_m1;
    1632          70 :     if (signe(n) < 0) x = ginv(x);
    1633          63 :     return Qp_sqrt(x);
    1634             :   }
    1635        1113 :   av = avma; p = gel(x,2);
    1636        1113 :   if (!signe(gel(x,4)))
    1637             :   {
    1638         203 :     if (signe(n) < 0) pari_err_INV("Qp_sqrtn", x);
    1639         203 :     q = divii(addis(n, valp(x)-1), n);
    1640         203 :     if (zetan) *zetan = gen_1;
    1641         203 :     avma = av; return zeropadic(p, itos(q));
    1642             :   }
    1643             :   /* treat the ramified part using logarithms */
    1644         910 :   e = Z_pvalrem(n, p, &q);
    1645         910 :   if (e) { x = Qp_sqrtn_ram(x,e); if (!x) return NULL; }
    1646         637 :   if (is_pm1(q))
    1647             :   { /* finished */
    1648          21 :     if (signe(q) < 0) x = ginv(x);
    1649          21 :     x = gerepileupto(av, x);
    1650          21 :     if (zetan)
    1651          28 :       *zetan = (e && absequaliu(p, 2))? gen_m1 /*-1 in Q_2*/
    1652          21 :                                    : gen_1;
    1653          21 :     return x;
    1654             :   }
    1655         616 :   tetpil = avma;
    1656             :   /* use hensel lift for unramified case */
    1657         616 :   x = Qp_sqrtn_unram(x, q, zetan);
    1658         616 :   if (!x) return NULL;
    1659         595 :   if (zetan)
    1660             :   {
    1661             :     GEN *gptr[2];
    1662          14 :     if (e && absequaliu(p, 2))/*-1 in Q_2*/
    1663             :     {
    1664           7 :       tetpil = avma; x = gcopy(x); *zetan = gneg(*zetan);
    1665             :     }
    1666          14 :     gptr[0] = &x; gptr[1] = zetan;
    1667          14 :     gerepilemanysp(av,tetpil,gptr,2);
    1668          14 :     return x;
    1669             :   }
    1670         581 :   return gerepile(av,tetpil,x);
    1671             : }
    1672             : 
    1673             : GEN
    1674         978 : sqrtnint(GEN a, long n)
    1675             : {
    1676         978 :   pari_sp ltop = avma;
    1677             :   GEN x, b, q;
    1678             :   long s, k, e;
    1679         978 :   const ulong nm1 = n - 1;
    1680         978 :   if (typ(a) != t_INT) pari_err_TYPE("sqrtnint",a);
    1681         978 :   if (n <= 0) pari_err_DOMAIN("sqrtnint", "n", "<=", gen_0, stoi(n));
    1682         971 :   if (n == 1) return icopy(a);
    1683         971 :   s = signe(a);
    1684         971 :   if (s < 0) pari_err_DOMAIN("sqrtnint", "x", "<", gen_0, a);
    1685         964 :   if (!s) return gen_0;
    1686         964 :   if (lgefint(a) == 3) return utoi(usqrtn(itou(a), n));
    1687         535 :   e = expi(a); k = e/(2*n);
    1688         535 :   if (k == 0)
    1689             :   {
    1690             :     long flag;
    1691         291 :     if (n > e) {avma = ltop; return gen_1;}
    1692         291 :     flag = cmpii(a, powuu(3, n)); avma = ltop;
    1693         291 :     return (flag < 0) ? gen_2: stoi(3);
    1694             :   }
    1695         244 :   if (e < n*(BITS_IN_LONG - 1))
    1696             :   {
    1697             :     ulong s, xs, qs;
    1698         162 :     s = 1 + e/n; xs = 1UL << s;
    1699         162 :     qs = itou(shifti(a, -nm1*s));
    1700        1491 :     while (qs < xs) {
    1701        1174 :       xs -= (xs - qs + nm1)/n;
    1702        1174 :       q = divii(a, powuu(xs, nm1));
    1703        1174 :       if (lgefint(q) > 3) break;
    1704        1167 :       qs = itou(q);
    1705             :     }
    1706         162 :     return utoi(xs);
    1707             :   }
    1708          82 :   b = addui(1, shifti(a, -n*k));
    1709          82 :   x = shifti(addui(1, sqrtnint(b, n)), k);
    1710          82 :   q = divii(a, powiu(x, nm1));
    1711         309 :   while (cmpii(q, x) < 0) /* a priori one iteration, no GC necessary */
    1712             :   {
    1713         145 :     x = subii(x, divis(addui(nm1, subii(x, q)), n));
    1714         145 :     q = divii(a, powiu(x, nm1));
    1715             :   }
    1716          82 :   return gerepileuptoleaf(ltop, x);
    1717             : }
    1718             : 
    1719             : ulong
    1720         429 : usqrtn(ulong a, ulong n)
    1721             : {
    1722             :   ulong x, s, q;
    1723         429 :   const ulong nm1 = n - 1;
    1724         429 :   if (!n) pari_err_DOMAIN("sqrtnint", "n", "=", gen_0, utoi(n));
    1725         429 :   if (n == 1 || a == 0) return a;
    1726         429 :   s = 1 + expu(a)/n; x = 1UL << s;
    1727         429 :   q = (nm1*s >= BITS_IN_LONG)? 0: a >> (nm1*s);
    1728        1694 :   while (q < x) {
    1729             :     ulong X;
    1730         836 :     x -= (x - q + nm1)/n;
    1731         836 :     X = upowuu(x, nm1);
    1732         836 :     q = X? a/X: 0;
    1733             :   }
    1734         429 :   return x;
    1735             : }
    1736             : 
    1737             : static ulong
    1738      207695 : cubic_prec_mask(long n)
    1739             : {
    1740      207695 :   long a = n, i;
    1741      207695 :   ulong mask = 0;
    1742     1163483 :   for(i = 1;; i++, mask *= 3)
    1743             :   {
    1744     1163483 :     long c = a%3;
    1745     1163483 :     if (c) mask += 3 - c;
    1746     1163483 :     a = (a+2)/3;
    1747     1371178 :     if (a==1) return mask + upowuu(3, i);
    1748      955788 :   }
    1749             : }
    1750             : 
    1751             : /* cubic Newton iteration, |a|^(1/n), assuming a != 0 */
    1752             : GEN
    1753      290872 : sqrtnr_abs(GEN a, long n)
    1754             : {
    1755             :   pari_sp av;
    1756             :   GEN x, b;
    1757             :   long eold, n1, n2, prec, B, v;
    1758             :   ulong mask;
    1759             : 
    1760      290872 :   if (n == 1) return mpabs(a);
    1761      290313 :   if (n == 2) return sqrtr_abs(a);
    1762             : 
    1763      284604 :   prec = realprec(a);
    1764      284604 :   B = prec2nbits(prec);
    1765      284604 :   n1 = n+1;
    1766      284604 :   n2 = 2*n; av = avma;
    1767      284604 :   v = expo(a) / n;
    1768      284604 :   if (v) a = shiftr(a, -n*v);
    1769             : 
    1770      284604 :   b = rtor(a, LOWDEFAULTPREC);
    1771      284604 :   x = mpexp(divru(logr_abs(b), n));
    1772      284604 :   if (prec == LOWDEFAULTPREC)
    1773             :   {
    1774       84808 :     if (v) shiftr_inplace(x, v);
    1775       84808 :     return gerepileuptoleaf(av, x);
    1776             :   }
    1777      199796 :   mask = cubic_prec_mask(B + BITS_IN_LONG-1);
    1778      199796 :   eold = 1;
    1779             :   for(;;)
    1780             :   { /* reach BITS_IN_LONG */
    1781      960888 :     long enew = eold * 3;
    1782      960888 :     enew -= mask % 3;
    1783      960888 :     if (enew > BITS_IN_LONG) break; /* back up one step */
    1784      761092 :     mask /= 3;
    1785      761092 :     eold = enew;
    1786      761092 :   }
    1787             :   for(;;)
    1788             :   {
    1789      360604 :     long pr, enew = eold * 3;
    1790             :     GEN y, z;
    1791      360604 :     enew -= mask % 3;
    1792      360604 :     mask /= 3;
    1793      360604 :     pr = nbits2prec(enew);
    1794      360604 :     b = rtor(a, pr); setsigne(b,1);
    1795      360604 :     x = rtor(x, pr);
    1796      360604 :     y = subrr(powru(x, n), b);
    1797      360604 :     z = divrr(y, addrr(mulur(n1, y), mulur(n2, b)));
    1798      360604 :     shiftr_inplace(z,1);
    1799      360604 :     x = mulrr(x, subsr(1,z));
    1800      360604 :     if (mask == 1)
    1801             :     {
    1802      199796 :       if (v) shiftr_inplace(x, v);
    1803      199796 :       return gerepileuptoleaf(av, gprec_wtrunc(x,prec));
    1804             :     }
    1805      160808 :     eold = enew;
    1806      160808 :   }
    1807             : }
    1808             : 
    1809             : static void
    1810       10990 : shiftc_inplace(GEN z, long d)
    1811             : {
    1812       10990 :   shiftr_inplace(gel(z,1), d);
    1813       10990 :   shiftr_inplace(gel(z,2), d);
    1814       10990 : }
    1815             : 
    1816             : /* exp(2*Pi*I/n), same iteration as sqrtnr_abs, different initial point */
    1817             : static GEN
    1818       61507 : sqrtnof1(ulong n, long prec)
    1819             : {
    1820             :   pari_sp av;
    1821             :   GEN x;
    1822             :   long eold, n1, n2, B;
    1823             :   ulong mask;
    1824             : 
    1825       61507 :   B = prec2nbits(prec);
    1826       61507 :   n1 = n+1;
    1827       61507 :   n2 = 2*n; av = avma;
    1828             : 
    1829       61507 :   x = expIr(divru(Pi2n(1, LOWDEFAULTPREC), n));
    1830       61507 :   if (prec == LOWDEFAULTPREC) return gerepileupto(av, x);
    1831        7899 :   mask = cubic_prec_mask(B + BITS_IN_LONG-1);
    1832        7899 :   eold = 1;
    1833             :   for(;;)
    1834             :   { /* reach BITS_IN_LONG */
    1835       38696 :     long enew = eold * 3;
    1836       38696 :     enew -= mask % 3;
    1837       38696 :     if (enew > BITS_IN_LONG) break; /* back up one step */
    1838       30797 :     mask /= 3;
    1839       30797 :     eold = enew;
    1840       30797 :   }
    1841             :   for(;;)
    1842             :   {
    1843       10990 :     long pr, enew = eold * 3;
    1844             :     GEN y, z;
    1845       10990 :     enew -= mask % 3;
    1846       10990 :     mask /= 3;
    1847       10990 :     pr = nbits2prec(enew);
    1848       10990 :     x = cxtofp(x, pr);
    1849       10990 :     y = gsub(gpowgs(x, n), gen_1);
    1850       10990 :     z = gdiv(y, gaddgs(gmulsg(n1, y), n2));
    1851       10990 :     shiftc_inplace(z,1);
    1852       10990 :     x = gmul(x, gsubsg(1, z));
    1853       10990 :     if (mask == 1) return gerepileupto(av, gprec_w(x,prec));
    1854        3091 :     eold = enew;
    1855        3091 :   }
    1856             : }
    1857             : 
    1858             : /* exp(2iPi/d) */
    1859             : GEN
    1860      124995 : rootsof1u_cx(ulong n, long prec)
    1861             : {
    1862      124995 :   switch(n)
    1863             :   {
    1864          21 :     case 1: return gen_1;
    1865        5758 :     case 2: return gen_m1;
    1866       15985 :     case 4: return gen_I();
    1867             :     case 3: case 6: case 12:
    1868             :     {
    1869        1587 :       pari_sp av = avma;
    1870        1587 :       GEN a = (n == 3)? mkfrac(gen_m1,gen_2): ghalf;
    1871        1587 :       GEN sq3 = sqrtr_abs(utor(3, prec));
    1872        1587 :       shiftr_inplace(sq3, -1);
    1873        1587 :       a = (n == 12)? mkcomplex(sq3, a): mkcomplex(a, sq3);
    1874        1587 :       return gerepilecopy(av, a);
    1875             :     }
    1876             :     case 8:
    1877             :     {
    1878       40137 :       pari_sp av = avma;
    1879       40137 :       GEN sq2 = sqrtr_abs(utor(2, prec));
    1880       40137 :       shiftr_inplace(sq2,-1);
    1881       40137 :       return gerepilecopy(av, mkcomplex(sq2, sq2));
    1882             :     }
    1883             :   }
    1884       61507 :   return sqrtnof1(n, prec);
    1885             : }
    1886             : /* exp(2iPi/d), assume d a t_INT */
    1887             : GEN
    1888        3500 : rootsof1_cx(GEN d, long prec)
    1889             : {
    1890        3500 :   if (lgefint(d) == 3) return rootsof1u_cx((ulong)d[2], prec);
    1891           0 :   return expIr(divri(Pi2n(1,prec), d));
    1892             : }
    1893             : 
    1894             : GEN
    1895        2696 : gsqrtn(GEN x, GEN n, GEN *zetan, long prec)
    1896             : {
    1897             :   long i, lx, tx;
    1898             :   pari_sp av;
    1899             :   GEN y, z;
    1900        2696 :   if (typ(n)!=t_INT) pari_err_TYPE("sqrtn",n);
    1901        2696 :   if (!signe(n)) pari_err_DOMAIN("sqrtn", "n", "=", gen_0, n);
    1902        2696 :   if (is_pm1(n))
    1903             :   {
    1904          14 :     if (zetan) *zetan = gen_1;
    1905          14 :     return (signe(n) > 0)? gcopy(x): ginv(x);
    1906             :   }
    1907        2682 :   if (zetan) *zetan = gen_0;
    1908        2682 :   tx = typ(x);
    1909        2682 :   if (is_matvec_t(tx))
    1910             :   {
    1911           7 :     y = cgetg_copy(x, &lx);
    1912           7 :     for (i=1; i<lx; i++) gel(y,i) = gsqrtn(gel(x,i),n,NULL,prec);
    1913           7 :     return y;
    1914             :   }
    1915        2675 :   av = avma;
    1916        2675 :   switch(tx)
    1917             :   {
    1918             :   case t_INTMOD:
    1919             :     {
    1920          56 :       GEN p = gel(x,1), s;
    1921          56 :       z = gen_0;
    1922          56 :       y = cgetg(3,t_INTMOD);  gel(y,1) = icopy(p);
    1923          56 :       if (zetan) { z = cgetg(3,t_INTMOD); gel(z,1) = gel(y,1); }
    1924          56 :       s = Fp_sqrtn(gel(x,2),n,p,zetan);
    1925          35 :       if (!s) {
    1926          21 :         if (zetan) {avma=av; return gen_0;}
    1927          14 :         if (!BPSW_psp(p)) pari_err_PRIME("sqrtn [modulus]",p);
    1928           7 :         pari_err_SQRTN("gsqrtn",x);
    1929             :       }
    1930          14 :       gel(y,2) = s;
    1931          14 :       if (zetan) { gel(z,2) = *zetan; *zetan = z; }
    1932          14 :       return y;
    1933             :     }
    1934             : 
    1935             :   case t_PADIC:
    1936          56 :     y = Qp_sqrtn(x,n,zetan);
    1937          49 :     if (!y) {
    1938           7 :       if (zetan) return gen_0;
    1939           7 :       pari_err_SQRTN("gsqrtn",x);
    1940             :     }
    1941          42 :     return y;
    1942             : 
    1943          84 :   case t_FFELT: return FF_sqrtn(x,n,zetan);
    1944             : 
    1945             :   case t_INT: case t_FRAC: case t_REAL: case t_COMPLEX:
    1946        2297 :     i = precision(x); if (i) prec = i;
    1947        2297 :     if (isint1(x))
    1948           7 :       y = real_1(prec);
    1949        2290 :     else if (gequal0(x))
    1950             :     {
    1951             :       long b;
    1952          21 :       if (signe(n) < 0) pari_err_INV("gsqrtn",x);
    1953          21 :       if (isinexactreal(x))
    1954          14 :         b = sdivsi(gexpo(x), n);
    1955             :       else
    1956           7 :         b = -prec2nbits(prec);
    1957          21 :       if (typ(x) == t_COMPLEX)
    1958             :       {
    1959           7 :         y = cgetg(3,t_COMPLEX);
    1960           7 :         gel(y,1) = gel(y,2) = real_0_bit(b);
    1961             :       }
    1962             :       else
    1963          14 :         y = real_0_bit(b);
    1964             :     }
    1965             :     else
    1966             :     {
    1967        2269 :       long nn = itos_or_0(n);
    1968        2269 :       if (tx == t_INT) { x = itor(x,prec); tx = t_REAL; }
    1969        2269 :       if (nn > 0 && tx == t_REAL && signe(x) > 0)
    1970         939 :         y = sqrtnr(x, nn);
    1971             :       else
    1972        1330 :         y = gexp(gdiv(glog(x,prec), n), prec);
    1973        2269 :       y = gerepileupto(av, y);
    1974             :     }
    1975        2297 :     if (zetan) *zetan = rootsof1_cx(n, prec);
    1976        2297 :     return y;
    1977             : 
    1978             :   case t_QUAD:
    1979           7 :     return gsqrtn(quadtofp(x, prec), n, zetan, prec);
    1980             : 
    1981             :   default:
    1982         175 :     av = avma; if (!(y = toser_i(x))) break;
    1983         175 :     return gerepileupto(av, ser_powfrac(y, ginv(n), prec));
    1984             :   }
    1985           0 :   pari_err_TYPE("sqrtn",x);
    1986             :   return NULL;/* LCOV_EXCL_LINE */
    1987             : }
    1988             : 
    1989             : /********************************************************************/
    1990             : /**                                                                **/
    1991             : /**                             EXP(X) - 1                         **/
    1992             : /**                                                                **/
    1993             : /********************************************************************/
    1994             : /* exp(|x|) - 1, assume x != 0.
    1995             :  * For efficiency, x should be reduced mod log(2): if so, we have a < 0 */
    1996             : GEN
    1997     6873900 : exp1r_abs(GEN x)
    1998             : {
    1999     6873900 :   long l = realprec(x), a = expo(x), b = prec2nbits(l), L, i, n, m, B;
    2000             :   GEN y, p2, X;
    2001             :   pari_sp av;
    2002             :   double d;
    2003             : 
    2004     6873900 :   if (b + a <= 0) return mpabs(x);
    2005             : 
    2006     6868562 :   y = cgetr(l); av = avma;
    2007     6868562 :   B = b/3 + BITS_IN_LONG + (BITS_IN_LONG*BITS_IN_LONG)/ b;
    2008     6868562 :   d = a/2.; m = (long)(d + sqrt(d*d + B)); /* >= 0 */
    2009     6868562 :   if (m < (-a) * 0.1) m = 0; /* not worth it */
    2010     6868562 :   L = l + nbits2extraprec(m);
    2011             :  /* Multiplication is quadratic in this range (l is small, otherwise we
    2012             :   * use logAGM + Newton). Set Y = 2^(-e-a) x, compute truncated series
    2013             :   * sum x^k/k!: this costs roughly
    2014             :   *    m b^2 + sum_{k <= n} (k e + BITS_IN_LONG)^2
    2015             :   * bit operations with |x| <  2^(1+a), |Y| < 2^(1-e) , m = e+a and b bits of
    2016             :   * accuracy needed, so
    2017             :   *    B := (b / 3 + BITS_IN_LONG + BITS_IN_LONG^2 / b) ~ m(m-a)
    2018             :   * we want b ~ 3 m (m-a) or m~b+a hence
    2019             :   *     m = min( a/2 + sqrt(a^2/4 + B),  b + a )
    2020             :   * NB: e ~ (b/3)^(1/2) as b -> oo
    2021             :   *
    2022             :   * Truncate the sum at k = n (>= 1), the remainder is
    2023             :   *   sum_{k >= n+1} Y^k / k! < Y^(n+1) / (n+1)! (1-Y) < Y^(n+1) / n!
    2024             :   * We want Y^(n+1) / n! <= Y 2^-b, hence -n log_2 |Y| + log_2 n! >= b
    2025             :   *   log n! ~ (n + 1/2) log(n+1) - (n+1) + log(2Pi)/2,
    2026             :   * error bounded by 1/6(n+1) <= 1/12. Finally, we want
    2027             :   * n (-1/log(2) -log_2 |Y| + log_2(n+1)) >= b  */
    2028     6868562 :   b += m;
    2029     6868562 :   d = m-dbllog2(x)-1/LOG2; /* ~ -log_2 Y - 1/log(2) */
    2030     6868562 :   n = (long)(b / d);
    2031     6868562 :   if (n > 1)
    2032     6861531 :     n = (long)(b / (d + log2((double)n+1))); /* log~constant in small ranges */
    2033     6868562 :   while (n*(d+log2((double)n+1)) < b) n++; /* expect few corrections */
    2034             : 
    2035     6868562 :   X = rtor(x,L); shiftr_inplace(X, -m); setsigne(X, 1);
    2036     6868562 :   if (n == 1) p2 = X;
    2037             :   else
    2038             :   {
    2039     6868562 :     long s = 0, l1 = nbits2prec((long)(d + n + 16));
    2040     6868562 :     GEN unr = real_1(L);
    2041             :     pari_sp av2;
    2042             : 
    2043     6868562 :     p2 = cgetr(L); av2 = avma;
    2044    97762312 :     for (i=n; i>=2; i--, avma = av2)
    2045             :     { /* compute X^(n-1)/n! + ... + X/2 + 1 */
    2046             :       GEN p1, p3;
    2047    90893750 :       setprec(X,l1); p3 = divru(X,i);
    2048    90893750 :       l1 += dvmdsBIL(s - expo(p3), &s); if (l1>L) l1=L;
    2049    90893750 :       setprec(unr,l1); p1 = addrr_sign(unr,1, i == n? p3: mulrr(p3,p2),1);
    2050    90893750 :       setprec(p2,l1); affrr(p1,p2); /* p2 <- 1 + (X/i)*p2 */
    2051             :     }
    2052     6868562 :     setprec(X,L); p2 = mulrr(X,p2);
    2053             :   }
    2054             : 
    2055    75932392 :   for (i=1; i<=m; i++)
    2056             :   {
    2057    69063830 :     if (realprec(p2) > L) setprec(p2,L);
    2058    69063830 :     p2 = mulrr(p2, addsr(2,p2));
    2059             :   }
    2060     6868562 :   affrr_fixlg(p2,y); avma = av; return y;
    2061             : }
    2062             : 
    2063             : GEN
    2064        8942 : mpexpm1(GEN x)
    2065             : {
    2066        8942 :   long sx = signe(x);
    2067             :   GEN y, z;
    2068             :   pari_sp av;
    2069        8942 :   if (!sx) return real_0_bit(expo(x));
    2070        8935 :   if (sx > 0) return exp1r_abs(x);
    2071             :   /* compute exp(x) * (1 - exp(-x)) */
    2072        4135 :   av = avma; y = exp1r_abs(x);
    2073        4135 :   z = addsr(1, y); setsigne(z, -1);
    2074        4135 :   return gerepileupto(av, divrr(y, z));
    2075             : }
    2076             : 
    2077             : static GEN serexp(GEN x, long prec);
    2078             : GEN
    2079       10524 : gexpm1(GEN x, long prec)
    2080             : {
    2081       10524 :   switch(typ(x))
    2082             :   {
    2083        4025 :     case t_REAL: return mpexpm1(x);
    2084        4994 :     case t_COMPLEX: return cxexpm1(x,prec);
    2085          14 :     case t_PADIC: return gsubgs(Qp_exp(x), 1);
    2086             :     default:
    2087             :     {
    2088        1491 :       pari_sp av = avma;
    2089             :       long ey;
    2090             :       GEN y;
    2091        1491 :       if (!(y = toser_i(x))) break;
    2092        1470 :       ey = valp(y);
    2093        1470 :       if (ey < 0) pari_err_DOMAIN("expm1","valuation", "<", gen_0, x);
    2094        1470 :       if (gequal0(y)) return gcopy(y);
    2095        1463 :       if (ey)
    2096         336 :         return gerepileupto(av, gsubgs(serexp(y,prec), 1));
    2097             :       else
    2098             :       {
    2099        1127 :         GEN e1 = gexpm1(gel(y,2), prec), e = gaddgs(e1,1);
    2100        1127 :         y = gmul(e, serexp(serchop0(y),prec));
    2101        1127 :         gel(y,2) = e1;
    2102        1127 :         return gerepilecopy(av, y);
    2103             :       }
    2104             :     }
    2105             :   }
    2106          21 :   return trans_eval("expm1",gexpm1,x,prec);
    2107             : }
    2108             : /********************************************************************/
    2109             : /**                                                                **/
    2110             : /**                             EXP(X)                             **/
    2111             : /**                                                                **/
    2112             : /********************************************************************/
    2113             : 
    2114             : /* centermod(x, log(2)), set *sh to the quotient */
    2115             : static GEN
    2116     6813433 : modlog2(GEN x, long *sh)
    2117             : {
    2118     6813433 :   double d = rtodbl(x);
    2119     6813433 :   long q = (long) ((fabs(d) + (LOG2/2))/LOG2);
    2120     6813433 :   if (d > LOG2 * LONG_MAX)
    2121           0 :     pari_err_OVERFLOW("expo()"); /* avoid overflow in  q */
    2122     6813433 :   if (d < 0) q = -q;
    2123     6813433 :   *sh = q;
    2124     6813433 :   if (q) {
    2125     6517318 :     long l = realprec(x) + 1;
    2126     6517318 :     x = subrr(rtor(x,l), mulsr(q, mplog2(l)));
    2127     6517318 :     if (!signe(x)) return NULL;
    2128             :   }
    2129     6813433 :   return x;
    2130             : }
    2131             : 
    2132             : static GEN
    2133     6812742 : mpexp_basecase(GEN x)
    2134             : {
    2135     6812742 :   pari_sp av = avma;
    2136     6812742 :   long sh, l = realprec(x);
    2137             :   GEN y, z;
    2138             : 
    2139     6812742 :   y = modlog2(x, &sh);
    2140     6812742 :   if (!y) { avma = av; return real2n(sh, l); }
    2141     6812742 :   z = addsr(1, exp1r_abs(y));
    2142     6812742 :   if (signe(y) < 0) z = invr(z);
    2143     6812742 :   if (sh) {
    2144     6516627 :     shiftr_inplace(z, sh);
    2145     6516627 :     if (realprec(z) > l) z = rtor(z, l); /* spurious precision increase */
    2146             :   }
    2147             : #ifdef DEBUG
    2148             : {
    2149             :   GEN t = mplog(z), u = divrr(subrr(x, t),x);
    2150             :   if (signe(u) && expo(u) > 5-prec2nbits(minss(l,realprec(t))))
    2151             :     pari_err_BUG("exp");
    2152             : }
    2153             : #endif
    2154     6812742 :   return gerepileuptoleaf(av, z); /* NOT affrr, precision often increases */
    2155             : }
    2156             : 
    2157             : GEN
    2158     7048247 : mpexp(GEN x)
    2159             : {
    2160     7048247 :   const long s = 6; /*Initial steps using basecase*/
    2161     7048247 :   long i, p, l = realprec(x), sh;
    2162             :   GEN a, t, z;
    2163             :   ulong mask;
    2164             : 
    2165     7048247 :   if (l <= maxss(EXPNEWTON_LIMIT, (1L<<s) + 2))
    2166             :   {
    2167     7047556 :     if (!signe(x)) return mpexp0(x);
    2168     6812051 :     return mpexp_basecase(x);
    2169             :   }
    2170         691 :   z = cgetr(l); /* room for result */
    2171         691 :   x = modlog2(x, &sh);
    2172         691 :   if (!x) { avma = (pari_sp)(z+lg(z)); return real2n(sh, l); }
    2173         691 :   constpi(l); /* precompute for later logr_abs() */
    2174         691 :   mask = quadratic_prec_mask(prec2nbits(l)+BITS_IN_LONG);
    2175         691 :   for(i=0, p=1; i<s+TWOPOTBITS_IN_LONG; i++) { p <<= 1; if (mask & 1) p-=1; mask >>= 1; }
    2176         691 :   a = mpexp_basecase(rtor(x, nbits2prec(p)));
    2177         691 :   x = addrs(x,1);
    2178         691 :   if (realprec(x) < l+EXTRAPRECWORD) x = rtor(x, l+EXTRAPRECWORD);
    2179         691 :   a = rtor(a, l+EXTRAPRECWORD); /*append 0s */
    2180         691 :   t = NULL;
    2181             :   for(;;)
    2182             :   {
    2183         728 :     p <<= 1; if (mask & 1) p--;
    2184         728 :     mask >>= 1;
    2185         728 :     setprec(x, nbits2prec(p));
    2186         728 :     setprec(a, nbits2prec(p));
    2187         728 :     t = mulrr(a, subrr(x, logr_abs(a))); /* a (x - log(a)) */
    2188         728 :     if (mask == 1) break;
    2189          37 :     affrr(t, a); avma = (pari_sp)a;
    2190          37 :   }
    2191         691 :   affrr(t,z);
    2192         691 :   if (sh) shiftr_inplace(z, sh);
    2193         691 :   avma = (pari_sp)z; return z;
    2194             : }
    2195             : 
    2196             : static long
    2197       19934 : Qp_exp_prec(GEN x)
    2198             : {
    2199       19934 :   long k, e = valp(x), n = e + precp(x);
    2200       19934 :   GEN p = gel(x,2);
    2201       19934 :   int is2 = absequaliu(p,2);
    2202       19934 :   if (e < 1 || (e == 1 && is2)) return -1;
    2203       19906 :   if (is2)
    2204             :   {
    2205        6207 :     n--; e--; k = n/e;
    2206        6207 :     if (n%e == 0) k--;
    2207             :   }
    2208             :   else
    2209             :   { /* e > 0, n > 0 */
    2210       13699 :     GEN r, t = subiu(p, 1);
    2211       13699 :     k = itos(dvmdii(subiu(muliu(t,n), 1), subiu(muliu(t,e), 1), &r));
    2212       13699 :     if (!signe(r)) k--;
    2213             :   }
    2214       19906 :   return k;
    2215             : }
    2216             : 
    2217             : static GEN
    2218       21334 : Qp_exp_safe(GEN x)
    2219             : {
    2220             :   long k;
    2221             :   pari_sp av;
    2222             :   GEN y;
    2223             : 
    2224       21334 :   if (gequal0(x)) return gaddgs(x,1);
    2225       19864 :   k = Qp_exp_prec(x);
    2226       19864 :   if (k < 0) return NULL;
    2227       19857 :   av = avma;
    2228       19857 :   for (y=gen_1; k; k--) y = gaddsg(1, gdivgs(gmul(y,x), k));
    2229       19857 :   return gerepileupto(av, y);
    2230             : }
    2231             : 
    2232             : GEN
    2233       20872 : Qp_exp(GEN x)
    2234             : {
    2235       20872 :   GEN y = Qp_exp_safe(x);
    2236       20872 :   if (!y) pari_err_DOMAIN("gexp(t_PADIC)","argument","",gen_0,x);
    2237       20865 :   return y;
    2238             : }
    2239             : 
    2240             : static GEN
    2241          49 : cos_p(GEN x)
    2242             : {
    2243             :   long k;
    2244             :   pari_sp av;
    2245             :   GEN x2, y;
    2246             : 
    2247          49 :   if (gequal0(x)) return gaddgs(x,1);
    2248          28 :   k = Qp_exp_prec(x);
    2249          28 :   if (k < 0) return NULL;
    2250          21 :   av = avma; x2 = gsqr(x);
    2251          21 :   if (k & 1) k--;
    2252         105 :   for (y=gen_1; k; k-=2)
    2253             :   {
    2254          84 :     GEN t = gdiv(gmul(y,x2), muluu(k, k-1));
    2255          84 :     y = gsubsg(1, t);
    2256             :   }
    2257          21 :   return gerepileupto(av, y);
    2258             : }
    2259             : static GEN
    2260          63 : sin_p(GEN x)
    2261             : {
    2262             :   long k;
    2263             :   pari_sp av;
    2264             :   GEN x2, y;
    2265             : 
    2266          63 :   if (gequal0(x)) return gcopy(x);
    2267          42 :   k = Qp_exp_prec(x);
    2268          42 :   if (k < 0) return NULL;
    2269          28 :   av = avma; x2 = gsqr(x);
    2270          28 :   if (k & 1) k--;
    2271         133 :   for (y=gen_1; k; k-=2)
    2272             :   {
    2273         105 :     GEN t = gdiv(gmul(y,x2), muluu(k, k+1));
    2274         105 :     y = gsubsg(1, t);
    2275             :   }
    2276          28 :   return gerepileupto(av, gmul(y, x));
    2277             : }
    2278             : 
    2279             : static GEN
    2280     1178228 : cxexp(GEN x, long prec)
    2281             : {
    2282     1178228 :   GEN r, p1, p2, y = cgetg(3,t_COMPLEX);
    2283     1178228 :   pari_sp av = avma, tetpil;
    2284             :   long l;
    2285     1178228 :   l = precision(x); if (l > prec) prec = l;
    2286     1178228 :   r = gexp(gel(x,1),prec);
    2287     1178228 :   if (gequal0(r)) { gel(y,1) = r; gel(y,2) = r; return y; }
    2288     1178228 :   gsincos(gel(x,2),&p2,&p1,prec);
    2289     1178228 :   tetpil = avma;
    2290     1178228 :   gel(y,1) = gmul(r,p1);
    2291     1178228 :   gel(y,2) = gmul(r,p2);
    2292     1178228 :   gerepilecoeffssp(av,tetpil,y+1,2);
    2293     1178228 :   return y;
    2294             : }
    2295             : 
    2296             : /* given a t_SER x^v s(x), with s(0) != 0, return x^v(s - s(0)), shallow */
    2297             : GEN
    2298       16646 : serchop0(GEN s)
    2299             : {
    2300       16646 :   long i, l = lg(s);
    2301             :   GEN y;
    2302       16646 :   if (l == 2) return s;
    2303       16646 :   if (l == 3 && isexactzero(gel(s,2))) return s;
    2304       16646 :   y = cgetg(l, t_SER); y[1] = s[1];
    2305       16646 :   gel(y,2) = gen_0; for (i=3; i <l; i++) gel(y,i) = gel(s,i);
    2306       16646 :   return normalize(y);
    2307             : }
    2308             : 
    2309             : static GEN
    2310       30240 : serexp(GEN x, long prec)
    2311             : {
    2312             :   pari_sp av;
    2313             :   long i,j,lx,ly,ex,mi;
    2314             :   GEN p1,y,xd,yd;
    2315             : 
    2316       30240 :   ex = valp(x);
    2317       30240 :   if (ex < 0) pari_err_DOMAIN("exp","valuation", "<", gen_0, x);
    2318       30233 :   if (gequal0(x)) return gaddsg(1,x);
    2319       22939 :   lx = lg(x);
    2320       22939 :   if (ex)
    2321             :   {
    2322       15869 :     ly = lx+ex; y = cgetg(ly,t_SER);
    2323       15869 :     mi = lx-1; while (mi>=3 && isrationalzero(gel(x,mi))) mi--;
    2324       15869 :     mi += ex-2;
    2325       15869 :     y[1] = evalsigne(1) | _evalvalp(0) | evalvarn(varn(x));
    2326             :     /* zd[i] = coefficient of X^i in z */
    2327       15869 :     xd = x+2-ex; yd = y+2; ly -= 2;
    2328       15869 :     gel(yd,0) = gen_1;
    2329       15869 :     for (i=1; i<ex; i++) gel(yd,i) = gen_0;
    2330       79723 :     for (   ; i<ly; i++)
    2331             :     {
    2332       63854 :       av = avma; p1 = gen_0;
    2333     2057125 :       for (j=ex; j<=minss(i,mi); j++)
    2334     1993271 :         p1 = gadd(p1, gmulgs(gmul(gel(xd,j),gel(yd,i-j)),j));
    2335       63854 :       gel(yd,i) = gerepileupto(av, gdivgs(p1,i));
    2336             :     }
    2337       15869 :     return y;
    2338             :   }
    2339        7070 :   av = avma;
    2340        7070 :   return gerepileupto(av, gmul(gexp(gel(x,2),prec), serexp(serchop0(x),prec)));
    2341             : }
    2342             : 
    2343             : GEN
    2344     5128636 : gexp(GEN x, long prec)
    2345             : {
    2346     5128636 :   switch(typ(x))
    2347             :   {
    2348     3748529 :     case t_REAL: return mpexp(x);
    2349     1178228 :     case t_COMPLEX: return cxexp(x,prec);
    2350          42 :     case t_PADIC: return Qp_exp(x);
    2351             :     default:
    2352             :     {
    2353      201837 :       pari_sp av = avma;
    2354             :       GEN y;
    2355      201837 :       if (!(y = toser_i(x))) break;
    2356       21707 :       return gerepileupto(av, serexp(y,prec));
    2357             :     }
    2358             :   }
    2359      180130 :   return trans_eval("exp",gexp,x,prec);
    2360             : }
    2361             : 
    2362             : /********************************************************************/
    2363             : /**                                                                **/
    2364             : /**                           AGM(X, Y)                            **/
    2365             : /**                                                                **/
    2366             : /********************************************************************/
    2367             : static int
    2368     3664038 : agmr_gap(GEN a, GEN b, long L)
    2369             : {
    2370     3664038 :   GEN d = subrr(b, a);
    2371     3664038 :   return (signe(d) && expo(d) - expo(b) >= L);
    2372             : }
    2373             : /* assume x > 0 */
    2374             : static GEN
    2375      266621 : agm1r_abs(GEN x)
    2376             : {
    2377      266621 :   long l = realprec(x), L = 5-prec2nbits(l);
    2378      266621 :   GEN a1, b1, y = cgetr(l);
    2379      266621 :   pari_sp av = avma;
    2380             : 
    2381      266621 :   a1 = addrr(real_1(l), x); shiftr_inplace(a1, -1);
    2382      266621 :   b1 = sqrtr_abs(x);
    2383     3930659 :   while (agmr_gap(a1,b1,L))
    2384             :   {
    2385     3397417 :     GEN a = a1;
    2386     3397417 :     a1 = addrr(a,b1); shiftr_inplace(a1, -1);
    2387     3397417 :     b1 = sqrtr_abs(mulrr(a,b1));
    2388             :   }
    2389      266621 :   affrr_fixlg(a1,y); avma = av; return y;
    2390             : }
    2391             : 
    2392             : struct agmcx_gap_t { long L, ex, cnt; };
    2393             : 
    2394             : static void
    2395       11267 : agmcx_init(GEN x, long *prec, struct agmcx_gap_t *S)
    2396             : {
    2397       11267 :   long l = precision(x);
    2398       11267 :   if (l) *prec = l;
    2399       11267 :   S->L = 1-prec2nbits(*prec);
    2400       11267 :   S->cnt = 0;
    2401       11267 :   S->ex = LONG_MAX;
    2402       11267 : }
    2403             : 
    2404             : static long
    2405       11267 : agmcx_a_b(GEN x, GEN *a1, GEN *b1, long prec)
    2406             : {
    2407       11267 :   long rotate = 0;
    2408       11267 :   if (gsigne(real_i(x))<0)
    2409             :   { /* Rotate by +/-Pi/2, so that the choice of the principal square
    2410             :      * root gives the optimal AGM. So a1 = +/-I*a1, b1=sqrt(-x). */
    2411         560 :     if (gsigne(imag_i(x))<0) { *a1=mulcxI(*a1);  rotate=-1; }
    2412         308 :     else                     { *a1=mulcxmI(*a1); rotate=1; }
    2413         560 :     x = gneg(x);
    2414             :   }
    2415       11267 :   *b1 = gsqrt(x, prec);
    2416       11267 :   return rotate;
    2417             : }
    2418             : /* return 0 if we must stop the AGM loop (a=b or a ~ b), 1 otherwise */
    2419             : static int
    2420      193742 : agmcx_gap(GEN a, GEN b, struct agmcx_gap_t *S)
    2421             : {
    2422      193742 :   GEN d = gsub(b, a);
    2423      193742 :   long ex = S->ex;
    2424      193742 :   S->ex = gexpo(d);
    2425      193742 :   if (gequal0(d) || S->ex - gexpo(b) < S->L) return 0;
    2426             :   /* if (S->ex >= ex) we're no longer making progress; twice in a row */
    2427      185659 :   if (S->ex < ex) S->cnt = 0;
    2428             :   else
    2429        6412 :     if (S->cnt++) return 0;
    2430      182475 :   return 1;
    2431             : }
    2432             : static GEN
    2433       11253 : agm1cx(GEN x, long prec)
    2434             : {
    2435             :   struct agmcx_gap_t S;
    2436             :   GEN a1, b1;
    2437       11253 :   pari_sp av = avma;
    2438             :   long rotate;
    2439       11253 :   agmcx_init(x, &prec, &S);
    2440       11253 :   a1 = gtofp(gmul2n(gadd(real_1(prec), x), -1), prec);
    2441       11253 :   rotate = agmcx_a_b(x, &a1, &b1, prec);
    2442      204909 :   while (agmcx_gap(a1,b1,&S))
    2443             :   {
    2444      182403 :     GEN a = a1;
    2445      182403 :     a1 = gmul2n(gadd(a,b1),-1);
    2446      182403 :     b1 = gsqrt(gmul(a,b1), prec);
    2447             :   }
    2448       11253 :   if (rotate) a1 = rotate>0 ? mulcxI(a1):mulcxmI(a1);
    2449       11253 :   return gerepilecopy(av,a1);
    2450             : }
    2451             : 
    2452             : GEN
    2453          14 : zellagmcx(GEN a0, GEN b0, GEN r, GEN t, long prec)
    2454             : {
    2455             :   struct agmcx_gap_t S;
    2456          14 :   pari_sp av = avma;
    2457          14 :   GEN x = gdiv(a0, b0), a1, b1;
    2458             :   long rotate;
    2459          14 :   agmcx_init(x, &prec, &S);
    2460          14 :   a1 = gtofp(gmul2n(gadd(real_1(prec), x), -1), prec);
    2461          14 :   r = gsqrt(gdiv(gmul(a1,gaddgs(r, 1)),gadd(r, x)), prec);
    2462          14 :   t = gmul(r, t);
    2463          14 :   rotate = agmcx_a_b(x, &a1, &b1, prec);
    2464         100 :   while (agmcx_gap(a1,b1,&S))
    2465             :   {
    2466          72 :     GEN a = a1, b = b1;
    2467          72 :     a1 = gmul2n(gadd(a,b),-1);
    2468          72 :     b1 = gsqrt(gmul(a,b), prec);
    2469          72 :     r = gsqrt(gdiv(gmul(a1,gaddgs(r, 1)),gadd(gmul(b, r), a )), prec);
    2470          72 :     t = gmul(r, t);
    2471             :   }
    2472          14 :   if (rotate) a1 = rotate>0 ? mulcxI(a1):mulcxmI(a1);
    2473          14 :   a1 = gmul(a1, b0);
    2474          14 :   t = gatan(gdiv(a1,t), prec);
    2475             :   /* send t to the fundamental domain if necessary */
    2476          14 :   if (gsigne(real_i(t))<0) t = gadd(t, mppi(prec));
    2477          14 :   return gerepileupto(av,gdiv(t,a1));
    2478             : }
    2479             : 
    2480             : static long
    2481          56 : ser_cmp_expo(GEN A, GEN B)
    2482             : {
    2483          56 :   long e = -(long)HIGHEXPOBIT, d = valp(B) - valp(A);
    2484          56 :   long i, la = lg(A), v = varn(B);
    2485       11256 :   for (i = 2; i < la; i++)
    2486             :   {
    2487       11200 :     GEN a = gel(A,i), b;
    2488             :     long ei;
    2489       11200 :     if (isexactzero(a)) continue;
    2490       11200 :     b = polcoeff_i(B, i-2 + d, v);
    2491       11200 :     ei = gexpo(a);
    2492       11200 :     if (!isexactzero(b)) ei -= gexpo(b);
    2493       11200 :     e = maxss(e, ei);
    2494             :   }
    2495          56 :   return e;
    2496             : }
    2497             : 
    2498             : static GEN
    2499          14 : ser_agm1(GEN y, long prec)
    2500             : {
    2501          14 :   GEN a1 = y, b1 = gen_1;
    2502          14 :   long l = lg(y)-2, l2 = 6-prec2nbits(prec), eold = LONG_MAX;
    2503             :   for(;;)
    2504             :   {
    2505          84 :     GEN a = a1, p1;
    2506          84 :     a1 = gmul2n(gadd(a,b1),-1);
    2507          84 :     b1 = gsqrt(gmul(a,b1), prec);
    2508          84 :     p1 = gsub(b1,a1);
    2509          84 :     if (isinexactreal(p1))
    2510             :     {
    2511          56 :       long e = ser_cmp_expo(p1, b1);
    2512          56 :       if (e < l2 || e >= eold) break;
    2513          49 :       eold = e;
    2514             :     }
    2515             :     else
    2516             :     {
    2517          28 :       long ep = valp(p1)-valp(b1);
    2518          28 :       if (ep >= l && gequal0(p1)) break;
    2519             :     }
    2520          70 :   }
    2521          14 :   return a1;
    2522             : }
    2523             : 
    2524             : /* agm(1,x) */
    2525             : static GEN
    2526        2676 : agm1(GEN x, long prec)
    2527             : {
    2528             :   GEN y;
    2529             :   pari_sp av;
    2530             : 
    2531        2676 :   if (gequal0(x)) return gcopy(x);
    2532        2676 :   switch(typ(x))
    2533             :   {
    2534             :     case t_INT:
    2535          28 :       if (!is_pm1(x)) break;
    2536          21 :       return (signe(x) > 0)? real_1(prec): real_0(prec);
    2537             : 
    2538        1381 :     case t_REAL: return signe(x) > 0? agm1r_abs(x): agm1cx(x, prec);
    2539             : 
    2540             :     case t_COMPLEX:
    2541        1134 :       if (gequal0(gel(x,2))) return agm1(gel(x,1), prec);
    2542        1134 :       return agm1cx(x, prec);
    2543             : 
    2544             :     case t_PADIC:
    2545             :     {
    2546          14 :       GEN a1 = x, b1 = gen_1;
    2547          14 :       long l = precp(x);
    2548          14 :       av = avma;
    2549             :       for(;;)
    2550             :       {
    2551          28 :         GEN a = a1, p1;
    2552             :         long ep;
    2553          28 :         a1 = gmul2n(gadd(a,b1),-1);
    2554          28 :         a = gmul(a,b1);
    2555          28 :         b1 = Qp_sqrt(a); if (!b1) pari_err_SQRTN("Qp_sqrt",a);
    2556          21 :         p1 = gsub(b1,a1); ep = valp(p1)-valp(b1);
    2557          21 :         if (ep<=0) { b1 = gneg_i(b1); p1 = gsub(b1,a1); ep=valp(p1)-valp(b1); }
    2558          21 :         if (ep >= l || gequal0(p1)) return gerepilecopy(av,a1);
    2559          14 :       }
    2560             :     }
    2561             : 
    2562             :     default:
    2563         119 :       av = avma; if (!(y = toser_i(x))) break;
    2564          14 :       return gerepilecopy(av, ser_agm1(y, prec));
    2565             :   }
    2566         112 :   return trans_eval("agm",agm1,x,prec);
    2567             : }
    2568             : 
    2569             : GEN
    2570        2564 : agm(GEN x, GEN y, long prec)
    2571             : {
    2572             :   pari_sp av;
    2573        2564 :   if (is_matvec_t(typ(y)))
    2574             :   {
    2575          14 :     if (is_matvec_t(typ(x))) pari_err_TYPE2("agm",x,y);
    2576           7 :     swap(x, y);
    2577             :   }
    2578        2557 :   if (gequal0(y)) return gcopy(y);
    2579        2557 :   av = avma;
    2580        2557 :   return gerepileupto(av, gmul(y, agm1(gdiv(x,y), prec)));
    2581             : }
    2582             : 
    2583             : /********************************************************************/
    2584             : /**                                                                **/
    2585             : /**                             LOG(X)                             **/
    2586             : /**                                                                **/
    2587             : /********************************************************************/
    2588             : /* atanh(u/v) using binary splitting */
    2589             : static GEN
    2590        4434 : atanhQ_split(ulong u, ulong v, long prec)
    2591             : {
    2592             :   long i, nmax;
    2593        4434 :   GEN u2 = sqru(u), v2 = sqru(v);
    2594        4434 :   double d = ((double)v) / u;
    2595             :   struct abpq_res R;
    2596             :   struct abpq A;
    2597             :   /* satisfies (2n+1) (v/u)^2n > 2^bitprec */
    2598        4434 :   nmax = (long)(prec2nbits(prec) / (2*log2(d)));
    2599        4434 :   abpq_init(&A, nmax);
    2600        4434 :   A.a[0] = A.b[0] = gen_1;
    2601        4434 :   A.p[0] = utoipos(u);
    2602        4434 :   A.q[0] = utoipos(v);
    2603      376390 :   for (i = 1; i <= nmax; i++)
    2604             :   {
    2605      371956 :     A.a[i] = gen_1;
    2606      371956 :     A.b[i] = utoipos((i<<1)+1);
    2607      371956 :     A.p[i] = u2;
    2608      371956 :     A.q[i] = v2;
    2609             :   }
    2610        4434 :   abpq_sum(&R, 0, nmax, &A);
    2611        4434 :   return rdivii(R.T, mulii(R.B,R.Q),prec);
    2612             : }
    2613             : /* log(2) = 10*atanh(1/17)+4*atanh(13/499); faster than logagmr_abs()
    2614             :  * and Pi/2M(1,4/2^n) ~ n log(2) */
    2615             : static GEN
    2616        2217 : log2_split(long prec)
    2617             : {
    2618        2217 :   GEN u = atanhQ_split(1, 17, prec);
    2619        2217 :   GEN v = atanhQ_split(13, 499, prec);
    2620        2217 :   shiftr_inplace(v, 2);
    2621        2217 :   return addrr(mulur(10, u), v);
    2622             : }
    2623             : GEN
    2624     9355437 : constlog2(long prec)
    2625             : {
    2626             :   pari_sp av;
    2627             :   GEN tmp;
    2628     9355437 :   if (glog2 && realprec(glog2) >= prec) return glog2;
    2629             : 
    2630        2217 :   tmp = cgetr_block(prec);
    2631        2217 :   av = avma;
    2632        2217 :   affrr(log2_split(prec+EXTRAPRECWORD), tmp);
    2633        2217 :   swap_clone(&glog2,tmp);
    2634        2217 :   avma = av; return glog2;
    2635             : }
    2636             : 
    2637             : GEN
    2638     9355437 : mplog2(long prec) { return rtor(constlog2(prec), prec); }
    2639             : 
    2640             : static GEN
    2641      295683 : logagmr_abs(GEN q)
    2642             : {
    2643      295683 :   long prec = realprec(q), lim, e = expo(q);
    2644             :   GEN z, y, Q, _4ovQ;
    2645             :   pari_sp av;
    2646             : 
    2647      295683 :   if (absrnz_equal2n(q)) return e? mulsr(e, mplog2(prec)): real_0(prec);
    2648      265282 :   z = cgetr(prec); av = avma; incrprec(prec);
    2649      265282 :   lim = prec2nbits(prec) >> 1;
    2650      265282 :   Q = rtor(q,prec);
    2651      265282 :   shiftr_inplace(Q,lim-e); setsigne(Q,1);
    2652             : 
    2653      265282 :   _4ovQ = invr(Q); shiftr_inplace(_4ovQ, 2); /* 4/Q */
    2654             :   /* Pi / 2agm(1, 4/Q) ~ log(Q), q = Q * 2^(e-lim) */
    2655      265282 :   y = divrr(Pi2n(-1, prec), agm1r_abs(_4ovQ));
    2656      265282 :   y = addrr(y, mulsr(e - lim, mplog2(prec)));
    2657      265282 :   affrr_fixlg(y, z); avma = av; return z;
    2658             : }
    2659             : /*return log(|x|), assuming x != 0 */
    2660             : GEN
    2661     2951863 : logr_abs(GEN X)
    2662             : {
    2663             :   pari_sp ltop;
    2664     2951863 :   long EX, L, m, k, a, b, l = realprec(X);
    2665             :   GEN z, x, y;
    2666             :   ulong u;
    2667             :   double d;
    2668             : 
    2669     2951863 :   if (l > LOGAGM_LIMIT) return logagmr_abs(X);
    2670             : 
    2671             :  /* Assuming 1 < x < 2, we want delta = x-1, 1-x/2, 1-1/x, or 2/x-1 small.
    2672             :   * We have 2/x-1 > 1-x/2, 1-1/x < x-1. So one should be choosing between
    2673             :   * 1-1/x and 1-x/2 ( crossover sqrt(2), worse ~ 0.29 ). To avoid an inverse,
    2674             :   * we choose between x-1 and 1-x/2 ( crossover 4/3, worse ~ 0.33 ) */
    2675     2656180 :   EX = expo(X);
    2676     2656180 :   u = uel(X,2);
    2677     2656180 :   k = 2;
    2678     2656180 :   if (u > (~0UL / 3) * 2) { /* choose 1-x/2 */
    2679     1574664 :     EX++; u = ~u;
    2680     1574664 :     while (!u && ++k < l) { u = uel(X,k); u = ~u; }
    2681             :   } else { /* choose x - 1 */
    2682     1081516 :     u &= ~HIGHBIT; /* u - HIGHBIT, assuming HIGHBIT set */
    2683     1081516 :     while (!u && ++k < l) u = uel(X,k);
    2684             :   }
    2685     2656180 :   if (k == l) return EX? mulsr(EX, mplog2(l)): real_0(l);
    2686     2561758 :   z = cgetr(EX? l: l - (k-2)); ltop = avma;
    2687             : 
    2688     2561758 :   a = prec2nbits(k) + bfffo(u); /* ~ -log2 |1-x| */
    2689             :  /* Multiplication is quadratic in this range (l is small, otherwise we
    2690             :   * use AGM). Set Y = x^(1/2^m), y = (Y - 1) / (Y + 1) and compute truncated
    2691             :   * series sum y^(2k+1)/(2k+1): the costs is less than
    2692             :   *    m b^2 + sum_{k <= n} ((2k+1) e + BITS_IN_LONG)^2
    2693             :   * bit operations with |x-1| <  2^(1-a), |Y| < 2^(1-e) , m = e-a and b bits of
    2694             :   * accuracy needed (+ BITS_IN_LONG since bit accuracies increase by
    2695             :   * increments of BITS_IN_LONG), so
    2696             :   * 4n^3/3 e^2 + n^2 2e BITS_IN_LONG+ n BITS_IN_LONG ~ m b^2, with n ~ b/2e
    2697             :   * or b/6e + BITS_IN_LONG/2e + BITS_IN_LONG/2be ~ m
    2698             :   *    B := (b / 6 + BITS_IN_LONG/2 + BITS_IN_LONG^2 / 2b) ~ m(m+a)
    2699             :   *     m = min( -a/2 + sqrt(a^2/4 + B),  b - a )
    2700             :   * NB: e ~ (b/6)^(1/2) as b -> oo */
    2701     2561758 :   L = l+1;
    2702     2561758 :   b = prec2nbits(L - (k-2)); /* take loss of accuracy into account */
    2703             :   /* instead of the above pessimistic estimate for the cost of the sum, use
    2704             :    * optimistic estimate (BITS_IN_LONG -> 0) */
    2705     2561758 :   d = -a/2.; m = (long)(d + sqrt(d*d + b/6)); /* >= 0 */
    2706             : 
    2707     2561758 :   if (m > b-a) m = b-a;
    2708     2561758 :   if (m < 0.2*a) m = 0; else L += nbits2extraprec(m);
    2709     2561758 :   x = rtor(X,L);
    2710     2561758 :   setsigne(x,1); shiftr_inplace(x,-EX);
    2711             :   /* 2/3 < x < 4/3 */
    2712     2561758 :   for (k=1; k<=m; k++) x = sqrtr_abs(x);
    2713             : 
    2714     2561758 :   y = divrr(subrs(x,1), addrs(x,1)); /* = (x-1) / (x+1), close to 0 */
    2715     2561758 :   L = realprec(y); /* should be ~ l+1 - (k-2) */
    2716             :   /* log(x) = log(1+y) - log(1-y) = 2 sum_{k odd} y^k / k
    2717             :    * Truncate the sum at k = 2n+1, the remainder is
    2718             :    *   2 sum_{k >= 2n+3} y^k / k < 2y^(2n+3) / (2n+3)(1-y) < y^(2n+3)
    2719             :    * We want y^(2n+3) < y 2^(-prec2nbits(L)), hence
    2720             :    *   n+1 > -prec2nbits(L) /-log_2(y^2) */
    2721     2561758 :   d = -2*dbllog2r(y); /* ~ -log_2(y^2) */
    2722     2561758 :   k = (long)(2*(prec2nbits(L) / d));
    2723     2561758 :   k |= 1;
    2724     2561758 :   if (k >= 3)
    2725             :   {
    2726     2558482 :     GEN S, T, y2 = sqrr(y), unr = real_1(L);
    2727     2558482 :     pari_sp av = avma;
    2728     2558482 :     long s = 0, incs = (long)d, l1 = nbits2prec((long)d);
    2729     2558482 :     S = x;
    2730     2558482 :     setprec(S,  l1);
    2731     2558482 :     setprec(unr,l1); affrr(divru(unr,k), S); /* destroy x, not needed anymore */
    2732    42183005 :     for (k -= 2;; k -= 2) /* k = 2n+1, ..., 1 */
    2733             :     { /* S = y^(2n+1-k)/(2n+1) + ... + 1 / k */
    2734    42183005 :       setprec(y2, l1); T = mulrr(S,y2);
    2735    42183005 :       if (k == 1) break;
    2736             : 
    2737    39624523 :       l1 += dvmdsBIL(s + incs, &s); if (l1>L) l1=L;
    2738    39624523 :       setprec(S, l1);
    2739    39624523 :       setprec(unr,l1);
    2740    39624523 :       affrr(addrr(divru(unr, k), T), S); avma = av;
    2741    39624523 :     }
    2742             :     /* k = 1 special-cased for eficiency */
    2743     2558482 :     y = mulrr(y, addsr(1,T)); /* = log(X)/2 */
    2744             :   }
    2745     2561758 :   shiftr_inplace(y, m + 1);
    2746     2561758 :   if (EX) y = addrr(y, mulsr(EX, mplog2(l+1)));
    2747     2561758 :   affrr_fixlg(y, z); avma = ltop; return z;
    2748             : }
    2749             : 
    2750             : /* assume Im(q) != 0 and precision(q) >= prec. Compute log(q) with accuracy
    2751             :  * prec [disregard input accuracy] */
    2752             : GEN
    2753       10077 : logagmcx(GEN q, long prec)
    2754             : {
    2755       10077 :   GEN z = cgetc(prec), y, Q, a, b;
    2756             :   long lim, e, ea, eb;
    2757       10077 :   pari_sp av = avma;
    2758       10077 :   int neg = 0;
    2759             : 
    2760       10077 :   incrprec(prec);
    2761       10077 :   if (gsigne(gel(q,1)) < 0) { q = gneg(q); neg = 1; }
    2762       10077 :   lim = prec2nbits(prec) >> 1;
    2763       10077 :   Q = gtofp(q, prec);
    2764       10077 :   a = gel(Q,1);
    2765       10077 :   b = gel(Q,2);
    2766       10077 :   if (gequal0(a)) {
    2767           0 :     affrr_fixlg(logr_abs(b), gel(z,1));
    2768           0 :     y = Pi2n(-1, prec);
    2769           0 :     if (signe(b) < 0) setsigne(y, -1);
    2770           0 :     affrr_fixlg(y, gel(z,2)); avma = av; return z;
    2771             :   }
    2772       10077 :   ea = expo(a);
    2773       10077 :   eb = expo(b);
    2774       10077 :   e = ea <= eb ? lim - eb : lim - ea;
    2775       10077 :   shiftr_inplace(a, e);
    2776       10077 :   shiftr_inplace(b, e);
    2777             : 
    2778             :   /* Pi / 2agm(1, 4/Q) ~ log(Q), q = Q * 2^e */
    2779       10077 :   y = gdiv(Pi2n(-1, prec), agm1cx( gdivsg(4, Q), prec ));
    2780       10077 :   a = gel(y,1);
    2781       10077 :   b = gel(y,2);
    2782       10077 :   a = addrr(a, mulsr(-e, mplog2(prec)));
    2783       10077 :   if (realprec(a) <= LOWDEFAULTPREC) a = real_0_bit(expo(a));
    2784       17297 :   if (neg) b = gsigne(b) <= 0? gadd(b, mppi(prec))
    2785        7220 :                              : gsub(b, mppi(prec));
    2786       10077 :   affrr_fixlg(a, gel(z,1));
    2787       10077 :   affrr_fixlg(b, gel(z,2)); avma = av; return z;
    2788             : }
    2789             : 
    2790             : GEN
    2791      180881 : mplog(GEN x)
    2792             : {
    2793      180881 :   if (signe(x)<=0) pari_err_DOMAIN("mplog", "argument", "<=", gen_0, x);
    2794      180881 :   return logr_abs(x);
    2795             : }
    2796             : 
    2797             : /* pe = p^e, p prime, 0 < x < pe a t_INT coprime to p. Return the (p-1)-th
    2798             :  * root of 1 in (Z/pe)^* congruent to x mod p, resp x mod 4 if p = 2.
    2799             :  * Simplified form of Zp_sqrtnlift: 1/(p-1) is trivial to compute */
    2800             : GEN
    2801        1316 : Zp_teichmuller(GEN x, GEN p, long e, GEN pe)
    2802             : {
    2803             :   GEN q, z, p1;
    2804             :   pari_sp av;
    2805             :   ulong mask;
    2806        1316 :   if (absequaliu(p,2)) return (mod4(x) & 2)? subiu(pe,1): gen_1;
    2807        1316 :   if (e == 1) return icopy(x);
    2808        1316 :   av = avma;
    2809        1316 :   p1 = subiu(p, 1);
    2810        1316 :   mask = quadratic_prec_mask(e);
    2811        1316 :   q = p; z = remii(x, p);
    2812        7714 :   while (mask > 1)
    2813             :   { /* Newton iteration solving z^{1 - p} = 1, z = x (mod p) */
    2814        5082 :     GEN w, t, qold = q;
    2815        5082 :     if (mask <= 3) /* last iteration */
    2816        1316 :       q = pe;
    2817             :     else
    2818             :     {
    2819        3766 :       q = sqri(q);
    2820        3766 :       if (mask & 1) q = diviiexact(q, p);
    2821             :     }
    2822        5082 :     mask >>= 1;
    2823             :     /* q <= qold^2 */
    2824        5082 :     if (lgefint(q) == 3)
    2825             :     {
    2826        4935 :       ulong Z = uel(z,2), Q = uel(q,2), P1 = uel(p1,2);
    2827        4935 :       ulong W = (Q-1) / P1; /* -1/(p-1) + O(qold) */
    2828        4935 :       ulong T = Fl_mul(W, Fl_powu(Z,P1,Q) - 1, Q);
    2829        4935 :       Z = Fl_mul(Z, 1 + T, Q);
    2830        4935 :       z = utoi(Z);
    2831             :     }
    2832             :     else
    2833             :     {
    2834         147 :       w = diviiexact(subiu(qold,1),p1); /* -1/(p-1) + O(qold) */
    2835         147 :       t = Fp_mul(w, subiu(Fp_pow(z,p1,q), 1), q);
    2836         147 :       z = Fp_mul(z, addui(1,t), q);
    2837             :     }
    2838             :   }
    2839        1316 :   return gerepileuptoint(av, z);
    2840             : }
    2841             : 
    2842             : GEN
    2843        1197 : teichmullerinit(long p, long n)
    2844             : {
    2845             :   GEN t, pn, g, v;
    2846             :   ulong gp, tp;
    2847             :   long a, m;
    2848             : 
    2849        1197 :   if (p == 2) return mkvec(gen_1);
    2850        1197 :   if (!uisprime(p)) pari_err_PRIME("teichmullerinit",utoipos(p));
    2851             : 
    2852        1197 :   m = p >> 1; /* (p-1)/2 */
    2853        1197 :   tp= gp= pgener_Fl(p); /* order (p-1), gp^m = -1 */
    2854        1197 :   pn = powuu(p, n);
    2855        1197 :   v = cgetg(p, t_VEC);
    2856        1197 :   t = g = Zp_teichmuller(utoipos(gp), utoipos(p), n, pn);
    2857        1197 :   gel(v, 1) = gen_1;
    2858        1197 :   gel(v, p-1) = subiu(pn,1);
    2859        2842 :   for (a = 1; a < m; a++)
    2860             :   {
    2861        1645 :     gel(v, tp) = t;
    2862        1645 :     gel(v, p - tp) = Fp_neg(t, pn); /* g^(m+a) = -g^a */
    2863        1645 :     if (a < m-1)
    2864             :     {
    2865         896 :       t = Fp_mul(t, g, pn); /* g^(a+1) */
    2866         896 :       tp = Fl_mul(tp, gp, p); /* t mod p  */
    2867             :     }
    2868             :   }
    2869        1197 :   return v;
    2870             : }
    2871             : 
    2872             : /* tab from teichmullerinit or NULL */
    2873             : GEN
    2874         238 : teichmuller(GEN x, GEN tab)
    2875             : {
    2876             :   GEN p, q, y, z;
    2877         238 :   long n, tx = typ(x);
    2878             : 
    2879         238 :   if (!tab)
    2880             :   {
    2881         126 :     if (tx == t_VEC && lg(x) == 3)
    2882             :     {
    2883           7 :       p = gel(x,1);
    2884           7 :       q = gel(x,2);
    2885           7 :       if (typ(p) == t_INT && typ(q) == t_INT)
    2886           7 :         return teichmullerinit(itos(p), itos(q));
    2887             :     }
    2888             :   }
    2889         112 :   else if (typ(tab) != t_VEC) pari_err_TYPE("teichmuller",tab);
    2890         231 :   if (tx!=t_PADIC) pari_err_TYPE("teichmuller",x);
    2891         231 :   z = gel(x,4);
    2892         231 :   if (!signe(z)) return gcopy(x);
    2893         231 :   p = gel(x,2);
    2894         231 :   q = gel(x,3);
    2895         231 :   n = precp(x);
    2896         231 :   y = cgetg(5,t_PADIC);
    2897         231 :   y[1] = evalprecp(n) | _evalvalp(0);
    2898         231 :   gel(y,2) = icopy(p);
    2899         231 :   gel(y,3) = icopy(q);
    2900         231 :   if (tab)
    2901             :   {
    2902         112 :     ulong pp = itou_or_0(p);
    2903         112 :     if (lg(tab) != (long)pp) pari_err_TYPE("teichmuller",tab);
    2904         112 :     z = gel(tab, umodiu(z, pp));
    2905         112 :     if (typ(z) != t_INT) pari_err_TYPE("teichmuller",tab);
    2906         112 :     z = remii(z, q);
    2907             :   }
    2908             :   else
    2909         119 :     z = Zp_teichmuller(z, p, n, q);
    2910         231 :   gel(y,4) = z;
    2911         231 :   return y;
    2912             : }
    2913             : GEN
    2914           0 : teich(GEN x) { return teichmuller(x, NULL); }
    2915             : 
    2916             : GEN
    2917     2332009 : glog(GEN x, long prec)
    2918             : {
    2919             :   pari_sp av, tetpil;
    2920             :   GEN y, p1;
    2921             :   long l;
    2922             : 
    2923     2332009 :   switch(typ(x))
    2924             :   {
    2925             :     case t_REAL:
    2926     1544377 :       if (signe(x) >= 0)
    2927             :       {
    2928     1284345 :         if (!signe(x)) pari_err_DOMAIN("log", "argument", "=", gen_0, x);
    2929     1284331 :         return logr_abs(x);
    2930             :       }
    2931      260032 :       retmkcomplex(logr_abs(x), mppi(realprec(x)));
    2932             : 
    2933             :     case t_FRAC:
    2934             :     {
    2935             :       GEN a, b;
    2936             :       long e1, e2;
    2937      119450 :       av = avma;
    2938      119450 :       a = gel(x,1);
    2939      119450 :       b = gel(x,2);
    2940      119450 :       e1 = expi(subii(a,b)); e2 = expi(b);
    2941      119450 :       if (e2 > e1) prec += nbits2nlong(e2 - e1);
    2942      119450 :       x = fractor(x, prec);
    2943      119450 :       return gerepileupto(av, glog(x, prec));
    2944             :     }
    2945             :     case t_COMPLEX:
    2946      433987 :       if (ismpzero(gel(x,2))) return glog(gel(x,1), prec);
    2947      418866 :       l = precision(x); if (l > prec) prec = l;
    2948      418866 :       if (ismpzero(gel(x,1)))
    2949             :       {
    2950        3332 :         GEN a = gel(x,2), b;
    2951        3332 :         av = avma; b = Pi2n(-1,prec);
    2952        3332 :         if (gsigne(a) < 0) { setsigne(b, -1); a = gabs(a,prec); }
    2953        3332 :         a = isint1(a) ? gen_0: glog(a,prec);
    2954        3332 :         return gerepilecopy(av, mkcomplex(a, b));
    2955             :       }
    2956      415534 :       if (prec >= LOGAGMCX_LIMIT) return logagmcx(x, prec);
    2957      405653 :       y = cgetg(3,t_COMPLEX);
    2958      405653 :       gel(y,2) = garg(x,prec);
    2959      405653 :       av = avma; p1 = glog(cxnorm(x),prec); tetpil = avma;
    2960      405653 :       gel(y,1) = gerepile(av,tetpil,gmul2n(p1,-1)); return y;
    2961             : 
    2962         357 :     case t_PADIC: return Qp_log(x);
    2963             :     default:
    2964      233838 :       av = avma; if (!(y = toser_i(x))) break;
    2965          49 :       if (!signe(y)) pari_err_DOMAIN("log", "argument", "=", gen_0, x);
    2966          49 :       if (valp(y)) pari_err_DOMAIN("log", "series valuation", "!=", gen_0, x);
    2967          42 :       p1 = integser(gdiv(derivser(y), y)); /* log(y)' = y'/y */
    2968          42 :       if (!gequal1(gel(y,2))) p1 = gadd(p1, glog(gel(y,2),prec));
    2969          42 :       return gerepileupto(av, p1);
    2970             :   }
    2971      233789 :   return trans_eval("log",glog,x,prec);
    2972             : }
    2973             : /********************************************************************/
    2974             : /**                                                                **/
    2975             : /**                        SINE, COSINE                            **/
    2976             : /**                                                                **/
    2977             : /********************************************************************/
    2978             : 
    2979             : /* Reduce x0 mod Pi/2 to x in [-Pi/4, Pi/4]. Return cos(x)-1 */
    2980             : static GEN
    2981     4449489 : mpcosm1(GEN x, long *ptmod8)
    2982             : {
    2983     4449489 :   long a = expo(x), l = realprec(x), b, L, i, n, m, B;
    2984             :   GEN y, p2, x2;
    2985             :   double d;
    2986             : 
    2987     4449489 :   n = 0;
    2988     4449489 :   if (a >= 0)
    2989             :   {
    2990             :     long p;
    2991             :     GEN q;
    2992     3357328 :     if (a > 30)
    2993             :     {
    2994         621 :       GEN z, pitemp = Pi2n(-2, nbits2prec(a + 32));
    2995         621 :       z = addrr(x,pitemp); /* = x + Pi/4 */
    2996         621 :       if (expo(z) >= bit_prec(z) + 3) pari_err_PREC("mpcosm1");
    2997         621 :       shiftr_inplace(pitemp, 1);
    2998         621 :       q = floorr( divrr(z,pitemp) ); /* round ( x / (Pi/2) ) */
    2999         621 :       p = l+EXTRAPRECWORD; x = rtor(x,p);
    3000             :     } else {
    3001     3356707 :       q = stoi((long)floor(rtodbl(x) / (M_PI/2) + 0.5));
    3002     3356707 :       p = l;
    3003             :     }
    3004     3357328 :     if (signe(q))
    3005             :     {
    3006     3357328 :       x = subrr(x, mulir(q, Pi2n(-1,p))); /* x mod Pi/2  */
    3007     3357328 :       a = expo(x);
    3008     3357328 :       if (!signe(x) && a >= 0) pari_err_PREC("mpcosm1");
    3009     3357328 :       n = mod4(q); if (n && signe(q) < 0) n = 4 - n;
    3010             :     }
    3011             :   }
    3012             :   /* a < 0 */
    3013     4449489 :   b = signe(x); *ptmod8 = (b < 0)? 4 + n: n;
    3014     4449489 :   if (!b) return real_0_bit(expo(x)*2 - 1);
    3015             : 
    3016     4259286 :   b = prec2nbits(l);
    3017     4259286 :   if (b + 2*a <= 0) {
    3018      177341 :     y = sqrr(x); shiftr_inplace(y, -1); setsigne(y, -1);
    3019      177341 :     return y;
    3020             :   }
    3021             : 
    3022     4081945 :   y = cgetr(l);
    3023     4081945 :   B = b/6 + BITS_IN_LONG + (BITS_IN_LONG*BITS_IN_LONG/2)/ b;
    3024     4081945 :   d = a/2.; m = (long)(d + sqrt(d*d + B)); /* >= 0 ,*/
    3025     4081945 :   if (m < (-a) * 0.1) m = 0; /* not worth it */
    3026     4081945 :   L = l + nbits2extraprec(m);
    3027             : 
    3028     4081945 :   b += m;
    3029     4081945 :   d = 2.0 * (m-dbllog2r(x)-1/LOG2); /* ~ 2( - log_2 Y - 1/log(2) ) */
    3030     4081945 :   n = (long)(b / d);
    3031     4081945 :   if (n > 1)
    3032     4054186 :     n = (long)(b / (d + log2((double)n+1))); /* log~constant in small ranges */
    3033     4081945 :   while (n*(d+log2((double)n+1)) < b) n++; /* expect few corrections */
    3034             : 
    3035             :  /* Multiplication is quadratic in this range (l is small, otherwise we
    3036             :   * use logAGM + Newton). Set Y = 2^(-e-a) x, compute truncated series
    3037             :   * sum Y^2k/(2k)!: this costs roughly
    3038             :   *   m b^2 + sum_{k <= n} (2k e + BITS_IN_LONG)^2
    3039             :   *   ~ floor(b/2e) b^2 / 3  + m b^2
    3040             :   * bit operations with |x| <  2^(1+a), |Y| < 2^(1-e) , m = e+a and b bits of
    3041             :   * accuracy needed, so
    3042             :   *    B := ( b / 6 + BITS_IN_LONG + BITS_IN_LONG^2 / 2b) ~ m(m-a)
    3043             :   * we want b ~ 6 m (m-a) or m~b+a hence
    3044             :   *     m = min( a/2 + sqrt(a^2/4 + b/6),  b/2 + a )
    3045             :   * NB1: e ~ (b/6)^(1/2) or b/2.
    3046             :   * NB2: We use b/4 instead of b/6 in the formula above: hand-optimized...
    3047             :   *
    3048             :   * Truncate the sum at k = n (>= 1), the remainder is
    3049             :   * < sum_{k >= n+1} Y^2k / 2k! < Y^(2n+2) / (2n+2)!(1-Y^2) < Y^(2n+2)/(2n+1)!
    3050             :   * We want ... <= Y^2 2^-b, hence -2n log_2 |Y| + log_2 (2n+1)! >= b
    3051             :   *   log n! ~ (n + 1/2) log(n+1) - (n+1) + log(2Pi)/2,
    3052             :   * error bounded by 1/6(n+1) <= 1/12. Finally, we want
    3053             :   * 2n (-1/log(2) - log_2 |Y| + log_2(2n+2)) >= b  */
    3054     4081945 :   x = rtor(x, L); shiftr_inplace(x, -m); setsigne(x, 1);
    3055     4081945 :   x2 = sqrr(x);
    3056     4081945 :   if (n == 1) { p2 = x2; shiftr_inplace(p2, -1); setsigne(p2, -1); } /*-Y^2/2*/
    3057             :   else
    3058             :   {
    3059     4081945 :     GEN unr = real_1(L);
    3060             :     pari_sp av;
    3061     4081945 :     long s = 0, l1 = nbits2prec((long)(d + n + 16));
    3062             : 
    3063     4081945 :     p2 = cgetr(L); av = avma;
    3064    25356009 :     for (i=n; i>=2; i--)
    3065             :     {
    3066             :       GEN p1;
    3067    21274064 :       setprec(x2,l1); p1 = divrunu(x2, 2*i-1);
    3068    21274064 :       l1 += dvmdsBIL(s - expo(p1), &s); if (l1>L) l1=L;
    3069    21274064 :       if (i != n) p1 = mulrr(p1,p2);
    3070    21274064 :       setprec(unr,l1); p1 = addrr_sign(unr,1, p1,-signe(p1));
    3071    21274064 :       setprec(p2,l1); affrr(p1,p2); avma = av;
    3072             :     }
    3073     4081945 :     shiftr_inplace(p2, -1); togglesign(p2); /* p2 := -p2/2 */
    3074     4081945 :     setprec(x2,L); p2 = mulrr(x2,p2);
    3075             :   }
    3076             :   /* Now p2 = sum {1<= i <=n} (-1)^i x^(2i) / (2i)! ~ cos(x) - 1 */
    3077    37625209 :   for (i=1; i<=m; i++)
    3078             :   { /* p2 = cos(x)-1 --> cos(2x)-1 */
    3079    33543264 :     p2 = mulrr(p2, addsr(2,p2));
    3080    33543264 :     shiftr_inplace(p2, 1);
    3081    33543264 :     if ((i & 31) == 0) p2 = gerepileuptoleaf((pari_sp)y, p2);
    3082             :   }
    3083     4081945 :   affrr_fixlg(p2,y); return y;
    3084             : }
    3085             : 
    3086             : /* sqrt (|1 - (1+x)^2|) = sqrt(|x*(x+2)|). Sends cos(x)-1 to |sin(x)| */
    3087             : static GEN
    3088     2880791 : mpaut(GEN x)
    3089             : {
    3090     2880791 :   pari_sp av = avma;
    3091     2880791 :   GEN t = mulrr(x, addsr(2,x)); /* != 0 */
    3092     2880791 :   if (!signe(t)) return real_0_bit(expo(t) >> 1);
    3093     2690595 :   return gerepileuptoleaf(av, sqrtr_abs(t));
    3094             : }
    3095             : 
    3096             : /********************************************************************/
    3097             : /**                            COSINE                              **/
    3098             : /********************************************************************/
    3099             : 
    3100             : GEN
    3101     2805479 : mpcos(GEN x)
    3102             : {
    3103             :   long mod8;
    3104             :   pari_sp av;
    3105             :   GEN y,p1;
    3106             : 
    3107     2805479 :   if (!signe(x)) {
    3108        1024 :     long l = nbits2prec(-expo(x));
    3109        1024 :     if (l < LOWDEFAULTPREC) l = LOWDEFAULTPREC;
    3110        1024 :     return real_1(l);
    3111             :   }
    3112             : 
    3113     2804455 :   av = avma; p1 = mpcosm1(x,&mod8);
    3114     2804455 :   switch(mod8)
    3115             :   {
    3116      860915 :     case 0: case 4: y = addsr(1,p1); break;
    3117      677368 :     case 1: case 7: y = mpaut(p1); togglesign(y); break;
    3118      663025 :     case 2: case 6: y = subsr(-1,p1); break;
    3119      603147 :     default:        y = mpaut(p1); break; /* case 3: case 5: */
    3120             :   }
    3121     2804455 :   return gerepileuptoleaf(av, y);
    3122             : }
    3123             : 
    3124             : /* convert INT or FRAC to REAL, which is later reduced mod 2Pi : avoid
    3125             :  * cancellation */
    3126             : static GEN
    3127       12362 : tofp_safe(GEN x, long prec)
    3128             : {
    3129       35931 :   return (typ(x) == t_INT || gexpo(x) > 0)? gadd(x, real_0(prec))
    3130       14287 :                                           : fractor(x, prec);
    3131             : }
    3132             : 
    3133             : GEN
    3134      215266 : gcos(GEN x, long prec)
    3135             : {
    3136             :   pari_sp av;
    3137             :   GEN r, u, v, y, u1, v1;
    3138             :   long i;
    3139             : 
    3140      215266 :   switch(typ(x))
    3141             :   {
    3142      214601 :     case t_REAL: return mpcos(x);
    3143             :     case t_COMPLEX:
    3144          42 :       if (isintzero(gel(x,1))) return gcosh(gel(x,2), prec);
    3145          28 :       i = precision(x); if (i) prec = i;
    3146          28 :       y = cgetc(prec); av = avma;
    3147          28 :       r = gexp(gel(x,2),prec);
    3148          28 :       v1 = gmul2n(addrr(invr(r),r), -1); /* = cos(I*Im(x)) */
    3149          28 :       u1 = subrr(v1, r); /* = - I*sin(I*Im(x)) */
    3150          28 :       gsincos(gel(x,1),&u,&v,prec);
    3151          28 :       affrr_fixlg(gmul(v1,v), gel(y,1));
    3152          28 :       affrr_fixlg(gmul(u1,u), gel(y,2)); avma = av; return y;
    3153             : 
    3154             :     case t_INT: case t_FRAC:
    3155         546 :       y = cgetr(prec); av = avma;
    3156         546 :       affrr_fixlg(mpcos(tofp_safe(x,prec)), y); avma = av; return y;
    3157             : 
    3158          49 :     case t_PADIC: y = cos_p(x);
    3159          49 :       if (!y) pari_err_DOMAIN("gcos(t_PADIC)","argument","",gen_0,x);
    3160          42 :       return y;
    3161             : 
    3162             :     default:
    3163          28 :       av = avma; if (!(y = toser_i(x))) break;
    3164          21 :       if (gequal0(y)) return gerepileupto(av, gaddsg(1,y));
    3165          21 :       if (valp(y) < 0)
    3166           7 :         pari_err_DOMAIN("cos","valuation", "<", gen_0, x);
    3167          14 :       gsincos(y,&u,&v,prec);
    3168          14 :       return gerepilecopy(av,v);
    3169             :   }
    3170           7 :   return trans_eval("cos",gcos,x,prec);
    3171             : }
    3172             : /********************************************************************/
    3173             : /**                             SINE                               **/
    3174             : /********************************************************************/
    3175             : 
    3176             : GEN
    3177      281551 : mpsin(GEN x)
    3178             : {
    3179             :   long mod8;
    3180             :   pari_sp av;
    3181             :   GEN y,p1;
    3182             : 
    3183      281551 :   if (!signe(x)) return real_0_bit(expo(x));
    3184             : 
    3185      279351 :   av = avma; p1 = mpcosm1(x,&mod8);
    3186      279351 :   switch(mod8)
    3187             :   {
    3188      206608 :     case 0: case 6: y=mpaut(p1); break;
    3189       22550 :     case 1: case 5: y=addsr(1,p1); break;
    3190       27985 :     case 2: case 4: y=mpaut(p1); togglesign(y); break;
    3191       22208 :     default:        y=subsr(-1,p1); break; /* case 3: case 7: */
    3192             :   }
    3193      279351 :   return gerepileuptoleaf(av, y);
    3194             : }
    3195             : 
    3196             : GEN
    3197      285310 : gsin(GEN x, long prec)
    3198             : {
    3199             :   pari_sp av;
    3200             :   GEN r, u, v, y, v1, u1;
    3201             :   long i;
    3202             : 
    3203      285310 :   switch(typ(x))
    3204             :   {
    3205      270197 :     case t_REAL: return mpsin(x);
    3206             :     case t_COMPLEX:
    3207        3612 :       if (isintzero(gel(x,1))) retmkcomplex(gen_0,gsinh(gel(x,2),prec));
    3208        3584 :       i = precision(x); if (i) prec = i;
    3209        3584 :       y = cgetc(prec); av = avma;
    3210        3584 :       r = gexp(gel(x,2),prec);
    3211        3584 :       v1 = gmul2n(addrr(invr(r),r), -1); /* = cos(I*Im(x)) */
    3212        3584 :       u1 = subrr(r, v1); /* = I*sin(I*Im(x)) */
    3213        3584 :       gsincos(gel(x,1),&u,&v,prec);
    3214        3584 :       affrr_fixlg(gmul(v1,u), gel(y,1));
    3215        3584 :       affrr_fixlg(gmul(u1,v), gel(y,2)); avma = av; return y;
    3216             : 
    3217             :     case t_INT: case t_FRAC:
    3218       11354 :       y = cgetr(prec); av = avma;
    3219       11354 :       affrr_fixlg(mpsin(tofp_safe(x,prec)), y); avma = av; return y;
    3220             : 
    3221          49 :     case t_PADIC: y = sin_p(x);
    3222          49 :       if (!y) pari_err_DOMAIN("gsin(t_PADIC)","argument","",gen_0,x);
    3223          42 :       return y;
    3224             : 
    3225             :     default:
    3226          98 :       av = avma; if (!(y = toser_i(x))) break;
    3227          91 :       if (gequal0(y)) return gerepilecopy(av, y);
    3228          91 :       if (valp(y) < 0)
    3229           7 :         pari_err_DOMAIN("sin","valuation", "<", gen_0, x);
    3230          84 :       gsincos(y,&u,&v,prec);
    3231          84 :       return gerepilecopy(av,u);
    3232             :   }
    3233           7 :   return trans_eval("sin",gsin,x,prec);
    3234             : }
    3235             : /********************************************************************/
    3236             : /**                       SINE, COSINE together                    **/
    3237             : /********************************************************************/
    3238             : 
    3239             : void
    3240     1383448 : mpsincos(GEN x, GEN *s, GEN *c)
    3241             : {
    3242             :   long mod8;
    3243             :   pari_sp av, tetpil;
    3244             :   GEN p1, *gptr[2];
    3245             : 
    3246     1383448 :   if (!signe(x))
    3247             :   {
    3248       22675 :     long e = expo(x);
    3249       22675 :     *s = real_0_bit(e);
    3250       22675 :     *c = e >= 0? real_0_bit(e): real_1_bit(-e);
    3251       45350 :     return;
    3252             :   }
    3253             : 
    3254     1360773 :   av=avma; p1=mpcosm1(x,&mod8); tetpil=avma;
    3255     1360773 :   switch(mod8)
    3256             :   {
    3257      308763 :     case 0: *c=addsr( 1,p1); *s=mpaut(p1); break;
    3258       88594 :     case 1: *s=addsr( 1,p1); *c=mpaut(p1); togglesign(*c); break;
    3259      289337 :     case 2: *c=subsr(-1,p1); *s=mpaut(p1); togglesign(*s); break;
    3260       98242 :     case 3: *s=subsr(-1,p1); *c=mpaut(p1); break;
    3261      223641 :     case 4: *c=addsr( 1,p1); *s=mpaut(p1); togglesign(*s); break;
    3262      107418 :     case 5: *s=addsr( 1,p1); *c=mpaut(p1); break;
    3263      145476 :     case 6: *c=subsr(-1,p1); *s=mpaut(p1); break;
    3264       99302 :     case 7: *s=subsr(-1,p1); *c=mpaut(p1); togglesign(*c); break;
    3265             :   }
    3266     1360773 :   gptr[0]=s; gptr[1]=c;
    3267     1360773 :   gerepilemanysp(av,tetpil,gptr,2);
    3268             : }
    3269             : 
    3270             : /* SINE and COSINE - 1 */
    3271             : void
    3272        4910 : mpsincosm1(GEN x, GEN *s, GEN *c)
    3273             : {
    3274             :   long mod8;
    3275             :   pari_sp av, tetpil;
    3276             :   GEN p1, *gptr[2];
    3277             : 
    3278        4910 :   if (!signe(x))
    3279             :   {
    3280           0 :     long e = expo(x);
    3281           0 :     *s = real_0_bit(e);
    3282           0 :     *c = real_0_bit(2*e-1);
    3283           0 :     return;
    3284             :   }
    3285        4910 :   av=avma; p1=mpcosm1(x,&mod8); tetpil=avma;
    3286        4910 :   switch(mod8)
    3287             :   {
    3288        4490 :     case 0: *c=rcopy(p1); *s=mpaut(p1); break;
    3289          42 :     case 1: *s=addsr(1,p1); *c=addrs(mpaut(p1),1); togglesign(*c); break;
    3290           0 :     case 2: *c=subsr(-2,p1); *s=mpaut(p1); togglesign(*s); break;
    3291           0 :     case 3: *s=subsr(-1,p1); *c=subrs(mpaut(p1),1); break;
    3292         287 :     case 4: *c=rcopy(p1); *s=mpaut(p1); togglesign(*s); break;
    3293          77 :     case 5: *s=addsr( 1,p1); *c=subrs(mpaut(p1),1); break;
    3294           7 :     case 6: *c=subsr(-2,p1); *s=mpaut(p1); break;
    3295           7 :     case 7: *s=subsr(-1,p1); *c=subsr(-1,mpaut(p1)); break;
    3296             :   }
    3297        4910 :   gptr[0]=s; gptr[1]=c;
    3298        4910 :   gerepilemanysp(av,tetpil,gptr,2);
    3299             : }
    3300             : 
    3301             : /* return exp(ix), x a t_REAL */
    3302             : GEN
    3303       89867 : expIr(GEN x)
    3304             : {
    3305       89867 :   pari_sp av = avma;
    3306       89867 :   GEN v = cgetg(3,t_COMPLEX);
    3307       89867 :   mpsincos(x, (GEN*)(v+2), (GEN*)(v+1));
    3308       89867 :   if (!signe(gel(v,2))) return gerepilecopy(av, gel(v,1));
    3309       87557 :   return v;
    3310             : }
    3311             : 
    3312             : /* return exp(ix)-1, x a t_REAL */
    3313             : static GEN
    3314        4910 : expm1_Ir(GEN x)
    3315             : {
    3316        4910 :   pari_sp av = avma;
    3317        4910 :   GEN v = cgetg(3,t_COMPLEX);
    3318        4910 :   mpsincosm1(x, (GEN*)(v+2), (GEN*)(v+1));
    3319        4910 :   if (!signe(gel(v,2))) return gerepilecopy(av, gel(v,1));
    3320        4910 :   return v;
    3321             : }
    3322             : 
    3323             : /* return exp(z)-1, z complex */
    3324             : GEN
    3325        5001 : cxexpm1(GEN z, long prec)
    3326             : {
    3327        5001 :   pari_sp av = avma;
    3328        5001 :   GEN X, Y, x = real_i(z), y = imag_i(z);
    3329        5001 :   long l = precision(z);
    3330        5001 :   if (l) prec = l;
    3331        5001 :   if (typ(x) != t_REAL) x = gtofp(x, prec);
    3332        5001 :   if (typ(y) != t_REAL) y = gtofp(y, prec);
    3333        5001 :   if (gequal0(y)) return mpexpm1(x);
    3334        4910 :   if (gequal0(x)) return expm1_Ir(y);
    3335        4826 :   X = mpexpm1(x); /* t_REAL */
    3336        4826 :   Y = expm1_Ir(y);
    3337             :   /* exp(x+iy) - 1 = (exp(x)-1)(exp(iy)-1) + exp(x)-1 + exp(iy)-1 */
    3338        4826 :   return gerepileupto(av, gadd(gadd(X,Y), gmul(X,Y)));
    3339             : }
    3340             : 
    3341             : void
    3342     1182792 : gsincos(GEN x, GEN *s, GEN *c, long prec)
    3343             : {
    3344             :   long i, j, ex, ex2, lx, ly, mi;
    3345             :   pari_sp av, tetpil;
    3346             :   GEN y, r, u, v, u1, v1, p1, p2, p3, p4, ps, pc;
    3347             :   GEN *gptr[4];
    3348             : 
    3349     1182792 :   switch(typ(x))
    3350             :   {
    3351             :     case t_INT: case t_FRAC:
    3352         441 :       *s = cgetr(prec);
    3353         441 :       *c = cgetr(prec); av = avma;
    3354         441 :       mpsincos(tofp_safe(x, prec), &ps, &pc);
    3355         441 :       affrr_fixlg(ps,*s);
    3356     1183233 :       affrr_fixlg(pc,*c); avma = av; return;
    3357             : 
    3358             :     case t_REAL:
    3359     1181959 :       mpsincos(x,s,c); return;
    3360             : 
    3361             :     case t_COMPLEX:
    3362         154 :       i = precision(x); if (i) prec = i;
    3363         154 :       ps = cgetc(prec); *s = ps;
    3364         154 :       pc = cgetc(prec); *c = pc; av = avma;
    3365         154 :       r = gexp(gel(x,2),prec);
    3366         154 :       v1 = gmul2n(addrr(invr(r),r), -1); /* = cos(I*Im(x)) */
    3367         154 :       u1 = subrr(r, v1); /* = I*sin(I*Im(x)) */
    3368         154 :       gsincos(gel(x,1), &u,&v, prec);
    3369         154 :       affrr_fixlg(mulrr(v1,u), gel(ps,1));
    3370         154 :       affrr_fixlg(mulrr(u1,v), gel(ps,2));
    3371         154 :       affrr_fixlg(mulrr(v1,v), gel(pc,1));
    3372         154 :       affrr_fixlg(mulrr(u1,u), gel(pc,2)); togglesign(gel(pc,2));
    3373         154 :       avma = av; return;
    3374             : 
    3375             :     case t_QUAD:
    3376           0 :       av = avma; gsincos(quadtofp(x, prec), s, c, prec);
    3377           0 :       gerepileall(av, 2, s, c); return;
    3378             : 
    3379             :     default:
    3380         238 :       av = avma; if (!(y = toser_i(x))) break;
    3381         238 :       if (gequal0(y)) { *s = gerepilecopy(av,y); *c = gaddsg(1,*s); return; }
    3382             : 
    3383         238 :       ex = valp(y); lx = lg(y); ex2 = 2*ex+2;
    3384         238 :       if (ex < 0) pari_err_DOMAIN("gsincos","valuation", "<", gen_0, x);
    3385         238 :       if (ex2 > lx)
    3386             :       {
    3387          28 :         *s = x == y? gcopy(y): gerepilecopy(av, y); av = avma;
    3388          28 :         *c = gerepileupto(av, gsubsg(1, gdivgs(gsqr(y),2)));
    3389          28 :         return;
    3390             :       }
    3391         210 :       if (!ex)
    3392             :       {
    3393          49 :         gsincos(serchop0(y),&u,&v,prec);
    3394          49 :         gsincos(gel(y,2),&u1,&v1,prec);
    3395          49 :         p1 = gmul(v1,v);
    3396          49 :         p2 = gmul(u1,u);
    3397          49 :         p3 = gmul(v1,u);
    3398          49 :         p4 = gmul(u1,v); tetpil = avma;
    3399          49 :         *c = gsub(p1,p2);
    3400          49 :         *s = gadd(p3,p4);
    3401          49 :         gptr[0]=s; gptr[1]=c;
    3402          49 :         gerepilemanysp(av,tetpil,gptr,2);
    3403          49 :         return;
    3404             :       }
    3405             : 
    3406         161 :       ly = lx+2*ex;
    3407         161 :       mi = lx-1; while (mi>=3 && isrationalzero(gel(y,mi))) mi--;
    3408         161 :       mi += ex-2;
    3409         161 :       pc = cgetg(ly,t_SER); *c = pc;
    3410         161 :       ps = cgetg(lx,t_SER); *s = ps;
    3411         161 :       pc[1] = evalsigne(1) | _evalvalp(0) | evalvarn(varn(y));
    3412         161 :       gel(pc,2) = gen_1; ps[1] = y[1];
    3413         161 :       for (i=2; i<ex+2; i++) gel(ps,i) = gcopy(gel(y,i));
    3414         161 :       for (i=3; i< ex2; i++) gel(pc,i) = gen_0;
    3415        2331 :       for (i=ex2; i<ly; i++)
    3416             :       {
    3417        2170 :         long ii = i-ex;
    3418        2170 :         av = avma; p1 = gen_0;
    3419        4340 :         for (j=ex; j<=minss(ii-2,mi); j++)
    3420        2170 :           p1 = gadd(p1, gmulgs(gmul(gel(y,j-ex+2),gel(ps,ii-j)),j));
    3421        2170 :         gel(pc,i) = gerepileupto(av, gdivgs(p1,2-i));
    3422        2170 :         if (ii < lx)
    3423             :         {
    3424        2002 :           av = avma; p1 = gen_0;
    3425        3836 :           for (j=ex; j<=minss(i-ex2,mi); j++)
    3426        1834 :             p1 = gadd(p1,gmulgs(gmul(gel(y,j-ex+2),gel(pc,i-j)),j));
    3427        2002 :           p1 = gdivgs(p1,i-2);
    3428        2002 :           gel(ps,ii) = gerepileupto(av, gadd(p1,gel(y,ii)));
    3429             :         }
    3430             :       }
    3431         161 :       return;
    3432             :   }
    3433           0 :   pari_err_TYPE("gsincos",x);
    3434             : }
    3435             : 
    3436             : /********************************************************************/
    3437             : /**                                                                **/
    3438             : /**                              SINC                              **/
    3439             : /**                                                                **/
    3440             : /********************************************************************/
    3441             : static GEN
    3442      107436 : mpsinc(GEN x)
    3443             : {
    3444      107436 :   pari_sp av = avma;
    3445             :   GEN s, c;
    3446             : 
    3447      107436 :   if (!signe(x)) {
    3448           0 :     long l = nbits2prec(-expo(x));
    3449           0 :     if (l < LOWDEFAULTPREC) l = LOWDEFAULTPREC;
    3450           0 :     return real_1(l);
    3451             :   }
    3452             : 
    3453      107436 :   mpsincos(x,&s,&c);
    3454      107436 :   return gerepileuptoleaf(av, divrr(s,x));
    3455             : }
    3456             : 
    3457             : GEN
    3458      107506 : gsinc(GEN x, long prec)
    3459             : {
    3460             :   pari_sp av;
    3461             :   GEN r, u, v, y, u1, v1;
    3462             :   long i;
    3463             : 
    3464      107506 :   switch(typ(x))
    3465             :   {
    3466      107429 :     case t_REAL: return mpsinc(x);
    3467             :     case t_COMPLEX:
    3468          28 :       if (isintzero(gel(x,1)))
    3469             :       {
    3470          14 :         av = avma;
    3471          14 :         return gerepileuptoleaf(av,gdiv(gsinh(gel(x,2),prec),gel(x,2)));
    3472             :       }
    3473          14 :       i = precision(x); if (i) prec = i;
    3474          14 :       y = cgetc(prec); av = avma;
    3475          14 :       r = gexp(gel(x,2),prec);
    3476          14 :       v1 = gmul2n(addrr(invr(r),r), -1); /* = cos(I*Im(x)) */
    3477          14 :       u1 = subrr(r, v1); /* = I*sin(I*Im(x)) */
    3478          14 :       gsincos(gel(x,1),&u,&v,prec);
    3479          14 :       affc_fixlg(gdiv(mkcomplex(gmul(v1,u), gmul(u1,v)), x), y);
    3480          14 :       avma = av; return y;
    3481             : 
    3482             :     case t_INT:
    3483           7 :       if (!signe(x)) return real_1(prec); /*fall through*/
    3484             :     case t_FRAC:
    3485           7 :       y = cgetr(prec); av = avma;
    3486           7 :       affrr_fixlg(mpsinc(tofp_safe(x,prec)), y); avma = av; return y;
    3487             : 
    3488             :     case t_PADIC:
    3489          21 :       if (gequal0(x)) return cvtop(gen_1, gel(x,2), valp(x));
    3490          14 :       av = avma; y = sin_p(x);
    3491          14 :       if (!y) pari_err_DOMAIN("gsinc(t_PADIC)","argument","",gen_0,x);
    3492           7 :       return gerepileuptoleaf(av,gdiv(y,x));
    3493             : 
    3494             :     default:
    3495          14 :       av = avma; if (!(y = toser_i(x))) break;
    3496          14 :       if (gequal0(y)) return gerepileupto(av, gaddsg(1,y));
    3497          14 :       if (valp(y) < 0)
    3498           7 :         pari_err_DOMAIN("sinc","valuation", "<", gen_0, x);
    3499           7 :       gsincos(y,&u,&v,prec);
    3500           7 :       return gerepilecopy(av,gdiv(u,y));
    3501             :   }
    3502           0 :   return trans_eval("sinc",gsinc,x,prec);
    3503             : }
    3504             : 
    3505             : /********************************************************************/
    3506             : /**                                                                **/
    3507             : /**                     TANGENT and COTANGENT                      **/
    3508             : /**                                                                **/
    3509             : /********************************************************************/
    3510             : static GEN
    3511          35 : mptan(GEN x)
    3512             : {
    3513          35 :   pari_sp av = avma;
    3514             :   GEN s, c;
    3515             : 
    3516          35 :   mpsincos(x,&s,&c);
    3517          35 :   if (!signe(c))
    3518           0 :     pari_err_DOMAIN("tan", "argument", "=", strtoGENstr("Pi/2 + kPi"),x);
    3519          35 :   return gerepileuptoleaf(av, divrr(s,c));
    3520             : }
    3521             : 
    3522             : GEN
    3523         105 : gtan(GEN x, long prec)
    3524             : {
    3525             :   pari_sp av;
    3526             :   GEN y, s, c;
    3527             : 
    3528         105 :   switch(typ(x))
    3529             :   {
    3530          28 :     case t_REAL: return mptan(x);
    3531             : 
    3532             :     case t_COMPLEX: {
    3533          28 :       if (isintzero(gel(x,1))) retmkcomplex(gen_0,gtanh(gel(x,2),prec));
    3534          14 :       av = avma; y = mulcxmI(gtanh(mulcxI(x), prec)); /* tan x = -I th(I x) */
    3535          14 :       gel(y,1) = gcopy(gel(y,1));
    3536          14 :       return gerepileupto(av, y);
    3537             :     }
    3538             :     case t_INT: case t_FRAC:
    3539           7 :       y = cgetr(prec); av = avma;
    3540           7 :       affrr_fixlg(mptan(tofp_safe(x,prec)), y); avma = av; return y;
    3541             : 
    3542             :     case t_PADIC:
    3543          14 :       av = avma;
    3544          14 :       return gerepileupto(av, gdiv(gsin(x,prec), gcos(x,prec)));
    3545             : 
    3546             :     default:
    3547          28 :       av = avma; if (!(y = toser_i(x))) break;
    3548          21 :       if (gequal0(y)) return gerepilecopy(av, y);
    3549          21 :       if (valp(y) < 0)
    3550           7 :         pari_err_DOMAIN("tan","valuation", "<", gen_0, x);
    3551          14 :       gsincos(y,&s,&c,prec);
    3552          14 :       return gerepileupto(av, gdiv(s,c));
    3553             :   }
    3554           7 :   return trans_eval("tan",gtan,x,prec);
    3555             : }
    3556             : 
    3557             : static GEN
    3558          35 : mpcotan(GEN x)
    3559             : {
    3560          35 :   pari_sp av=avma, tetpil;
    3561             :   GEN s,c;
    3562             : 
    3563          35 :   mpsincos(x,&s,&c); tetpil=avma;
    3564          35 :   return gerepile(av,tetpil,divrr(c,s));
    3565             : }
    3566             : 
    3567             : GEN
    3568         154 : gcotan(GEN x, long prec)
    3569             : {
    3570             :   pari_sp av;
    3571             :   GEN y, s, c;
    3572             : 
    3573         154 :   switch(typ(x))
    3574             :   {
    3575             :     case t_REAL:
    3576          28 :       return mpcotan(x);
    3577             : 
    3578             :     case t_COMPLEX:
    3579          42 :       if (isintzero(gel(x,1))) {
    3580          21 :         GEN z = cgetg(3, t_COMPLEX);
    3581          21 :         gel(z,1) = gen_0;
    3582          21 :         av = avma;
    3583          21 :         gel(z,2) = gerepileupto(av, gneg(ginv(gtanh(gel(x,2),prec))));
    3584          21 :         return z;
    3585             :       }
    3586          21 :       av = avma;
    3587          21 :       gsincos(x,&s,&c,prec);
    3588          21 :       return gerepileupto(av, gdiv(c,s));
    3589             : 
    3590             :     case t_INT: case t_FRAC:
    3591           7 :       y = cgetr(prec); av = avma;
    3592           7 :       affrr_fixlg(mpcotan(tofp_safe(x,prec)), y); avma = av; return y;
    3593             : 
    3594             :     case t_PADIC:
    3595          14 :       av = avma;
    3596          14 :       return gerepileupto(av, gdiv(gcos(x,prec), gsin(x,prec)));
    3597             : 
    3598             :     default:
    3599          63 :       av = avma; if (!(y = toser_i(x))) break;
    3600          56 :       if (gequal0(y)) pari_err_DOMAIN("cotan", "argument", "=", gen_0, y);
    3601          56 :       if (valp(y) < 0) pari_err_DOMAIN("cotan","valuation", "<", gen_0, x);
    3602          49 :       gsincos(y,&s,&c,prec);
    3603          49 :       return gerepileupto(av, gdiv(c,s));
    3604             :   }
    3605           7 :   return trans_eval("cotan",gcotan,x,prec);
    3606             : }

Generated by: LCOV version 1.11