Code coverage tests

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

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

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

LCOV - code coverage report
Current view: top level - basemath - trans1.c (source / functions) Hit Total Coverage
Test: PARI/GP v2.12.1 lcov report (development 25819-e703fe1174) Lines: 2222 2288 97.1 %
Date: 2020-09-18 06:10:04 Functions: 163 165 98.8 %
Legend: Lines: hit not hit

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

Generated by: LCOV version 1.13