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.16.2 lcov report (development 29464-65cd88daf0) Lines: 2260 2325 97.2 %
Date: 2024-07-23 09:03:50 Functions: 166 167 99.4 %
Legend: Lines: hit not hit

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

Generated by: LCOV version 1.16