"Fossies" - the Fresh Open Source Software Archive

Member "pari-2.13.1/src/basemath/trans1.c" (14 Jan 2021, 105716 Bytes) of package /linux/misc/pari-2.13.1.tar.gz:


As a special service "Fossies" has tried to format the requested source page into HTML format using (guessed) C and C++ source code syntax highlighting (style: standard) with prefixed line numbers and code folding option. Alternatively you can here view or download the uninterpreted source code file. For more information about "trans1.c" see the Fossies "Dox" file reference documentation and the latest Fossies "Diffs" side-by-side code changes report: 2.13.0_vs_2.13.1.

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