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-bordeaux1.fr machine (x86_64 architecture), and agregate them in the final report:

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

LCOV - code coverage report
Current view: top level - basemath - trans1.c (source / functions) Hit Total Coverage
Test: PARI/GP v2.8.0 lcov report (development 17238-ce76349) Lines: 1687 1734 97.3 %
Date: 2014-12-20 Functions: 126 128 98.4 %
Legend: Lines: hit not hit | Branches: + taken - not taken # not executed Branches: 1054 1208 87.3 %

           Branch data     Line data    Source code
       1                 :            : /* Copyright (C) 2000  The PARI group.
       2                 :            : 
       3                 :            : This file is part of the PARI/GP package.
       4                 :            : 
       5                 :            : PARI/GP is free software; you can redistribute it and/or modify it under the
       6                 :            : terms of the GNU General Public License as published by the Free Software
       7                 :            : Foundation. It is distributed in the hope that it will be useful, but WITHOUT
       8                 :            : ANY WARRANTY WHATSOEVER.
       9                 :            : 
      10                 :            : Check the License for details. You should have received a copy of it, along
      11                 :            : with the package; see the file 'COPYING'. If not, write to the Free Software
      12                 :            : Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. */
      13                 :            : 
      14                 :            : /********************************************************************/
      15                 :            : /**                                                                **/
      16                 :            : /**                   TRANSCENDENTAL FUNCTIONS                     **/
      17                 :            : /**                                                                **/
      18                 :            : /********************************************************************/
      19                 :            : #include "pari.h"
      20                 :            : #include "paripriv.h"
      21                 :            : 
      22                 :            : #ifdef LONG_IS_64BIT
      23                 :            : static const long SQRTVERYBIGINT = 3037000500L; /* ceil(sqrt(LONG_MAX)) */
      24                 :            : #else
      25                 :            : static const long SQRTVERYBIGINT = 46341L;
      26                 :            : #endif
      27                 :            : 
      28                 :            : static THREAD GEN gcatalan, geuler, glog2, gpi;
      29                 :            : void
      30                 :       1243 : pari_init_floats(void)
      31                 :            : {
      32                 :       1243 :   gcatalan = geuler = gpi = bernzone = glog2 = NULL;
      33                 :       1243 : }
      34                 :            : 
      35                 :            : void
      36                 :       1231 : pari_close_floats(void)
      37                 :            : {
      38         [ +  + ]:       1231 :   if (gcatalan) gunclone(gcatalan);
      39         [ +  + ]:       1231 :   if (geuler) gunclone(geuler);
      40         [ +  + ]:       1231 :   if (gpi) gunclone(gpi);
      41         [ +  + ]:       1231 :   if (bernzone) gunclone(bernzone);
      42         [ +  + ]:       1231 :   if (glog2) gunclone(glog2);
      43                 :       1231 : }
      44                 :            : 
      45                 :            : /********************************************************************/
      46                 :            : /**                   GENERIC BINARY SPLITTING                     **/
      47                 :            : /**                    (Haible, Papanikolaou)                      **/
      48                 :            : /********************************************************************/
      49                 :            : void
      50                 :       5443 : abpq_init(struct abpq *A, long n)
      51                 :            : {
      52                 :       5443 :   A->a = (GEN*)new_chunk(n+1);
      53                 :       5443 :   A->b = (GEN*)new_chunk(n+1);
      54                 :       5443 :   A->p = (GEN*)new_chunk(n+1);
      55                 :       5443 :   A->q = (GEN*)new_chunk(n+1);
      56                 :       5443 : }
      57                 :            : static GEN
      58                 :     454272 : mulii3(GEN a, GEN b, GEN c) { return mulii(mulii(a,b),c); }
      59                 :            : static GEN
      60                 :     134697 : mulii4(GEN a, GEN b, GEN c, GEN d) { return mulii(mulii(a,b),mulii(c,d)); }
      61                 :            : 
      62                 :            : /* T_{n1,n1+1}, given P = p[n1]p[n1+1] */
      63                 :            : static GEN
      64                 :     134697 : T2(struct abpq *A, long n1, GEN P)
      65                 :            : {
      66                 :     134697 :   GEN u1 = mulii4(A->a[n1], A->p[n1], A->b[n1+1], A->q[n1+1]);
      67                 :     134697 :   GEN u2 = mulii3(A->b[n1],A->a[n1+1], P);
      68                 :     134697 :   return addii(u1, u2);
      69                 :            : }
      70                 :            : 
      71                 :            : /* assume n2 > n1. Compute sum_{n1 <= n < n2} a/b(n) p/q(n1)... p/q(n) */
      72                 :            : void
      73                 :     264015 : abpq_sum(struct abpq_res *r, long n1, long n2, struct abpq *A)
      74                 :            : {
      75                 :            :   struct abpq_res L, R;
      76                 :            :   GEN u1, u2;
      77                 :            :   pari_sp av;
      78                 :            :   long n;
      79   [ +  +  +  + ]:     264015 :   switch(n2 - n1)
      80                 :            :   {
      81                 :            :     GEN b, p, q;
      82                 :            :     case 1:
      83                 :         32 :       r->P = A->p[n1];
      84                 :         32 :       r->Q = A->q[n1];
      85                 :         32 :       r->B = A->b[n1];
      86                 :         32 :       r->T = mulii(A->a[n1], A->p[n1]);
      87                 :            :       return;
      88                 :            :     case 2:
      89                 :      73694 :       r->P = mulii(A->p[n1], A->p[n1+1]);
      90                 :      73694 :       r->Q = mulii(A->q[n1], A->q[n1+1]);
      91                 :      73694 :       r->B = mulii(A->b[n1], A->b[n1+1]);
      92                 :      73694 :       av = avma;
      93                 :      73694 :       r->T = gerepileuptoint(av, T2(A, n1, r->P));
      94                 :            :       return;
      95                 :            : 
      96                 :            :     case 3:
      97                 :      61003 :       p = mulii(A->p[n1+1], A->p[n1+2]);
      98                 :      61003 :       q = mulii(A->q[n1+1], A->q[n1+2]);
      99                 :      61003 :       b = mulii(A->b[n1+1], A->b[n1+2]);
     100                 :      61003 :       r->P = mulii(A->p[n1], p);
     101                 :      61003 :       r->Q = mulii(A->q[n1], q);
     102                 :      61003 :       r->B = mulii(A->b[n1], b);
     103                 :      61003 :       av = avma;
     104                 :      61003 :       u1 = mulii3(b, q, A->a[n1]);
     105                 :      61003 :       u2 = mulii(A->b[n1], T2(A, n1+1, p));
     106                 :      61003 :       r->T = gerepileuptoint(av, mulii(A->p[n1], addii(u1, u2)));
     107                 :            :       return;
     108                 :            :   }
     109                 :            : 
     110                 :     129286 :   av = avma;
     111                 :     129286 :   n = (n1 + n2) >> 1;
     112                 :     129286 :   abpq_sum(&L, n1, n, A);
     113                 :     129286 :   abpq_sum(&R, n, n2, A);
     114                 :            : 
     115                 :     129286 :   r->P = mulii(L.P, R.P);
     116                 :     129286 :   r->Q = mulii(L.Q, R.Q);
     117                 :     129286 :   r->B = mulii(L.B, R.B);
     118                 :     129286 :   u1 = mulii3(R.B,R.Q,L.T);
     119                 :     129286 :   u2 = mulii3(L.B,L.P,R.T);
     120                 :     129286 :   r->T = addii(u1,u2);
     121                 :     129286 :   avma = av;
     122                 :     129286 :   r->P = icopy(r->P);
     123                 :     129286 :   r->Q = icopy(r->Q);
     124                 :     129286 :   r->B = icopy(r->B);
     125                 :     264015 :   r->T = icopy(r->T);
     126                 :            : }
     127                 :            : 
     128                 :            : /********************************************************************/
     129                 :            : /**                                                                **/
     130                 :            : /**                               PI                               **/
     131                 :            : /**                                                                **/
     132                 :            : /********************************************************************/
     133                 :            : /* replace *old clone by c. Protect against SIGINT */
     134                 :            : static void
     135                 :       3815 : swap_clone(GEN *old, GEN c)
     136                 :            : {
     137                 :       3815 :   GEN tmp = *old;
     138                 :       3815 :   *old = c;
     139         [ +  + ]:       3815 :   if (tmp) gunclone(tmp);
     140                 :       3815 : }
     141                 :            : 
     142                 :            : /*                         ----
     143                 :            :  *  53360 (640320)^(1/2)   \    (6n)! (545140134 n + 13591409)
     144                 :            :  *  -------------------- = /    ------------------------------
     145                 :            :  *        Pi               ----   (n!)^3 (3n)! (-640320)^(3n)
     146                 :            :  *                         n>=0
     147                 :            :  *
     148                 :            :  * Ramanujan's formula + binary splitting */
     149                 :            : static GEN
     150                 :       1690 : pi_ramanujan(long prec)
     151                 :            : {
     152                 :       1690 :   const ulong B = 545140134, A = 13591409, C = 640320;
     153                 :       1690 :   const double alpha2 = 47.11041314; /* 3log(C/12) / log(2) */
     154                 :            :   long n, nmax, prec2;
     155                 :            :   struct abpq_res R;
     156                 :            :   struct abpq S;
     157                 :            :   GEN D, u;
     158                 :            : 
     159                 :       1690 :   nmax = (long)(1 + prec2nbits(prec)/alpha2);
     160                 :            : #ifdef LONG_IS_64BIT
     161                 :       1446 :   D = utoipos(10939058860032000); /* C^3/24 */
     162                 :            : #else
     163                 :        244 :   D = uutoi(2546948,495419392);
     164                 :            : #endif
     165                 :       1690 :   abpq_init(&S, nmax);
     166                 :       1690 :   S.a[0] = utoipos(A);
     167                 :       1690 :   S.b[0] = S.p[0] = S.q[0] = gen_1;
     168         [ +  + ]:      27049 :   for (n = 1; n <= nmax; n++)
     169                 :            :   {
     170                 :      25359 :     S.a[n] = addiu(muluu(B, n), A);
     171                 :      25359 :     S.b[n] = gen_1;
     172                 :      25359 :     S.p[n] = mulis(muluu(6*n-5, 2*n-1), 1-6*n);
     173                 :      25359 :     S.q[n] = mulii(sqru(n), muliu(D,n));
     174                 :            :   }
     175                 :       1690 :   abpq_sum(&R, 0, nmax, &S); prec2 = prec+EXTRAPRECWORD;
     176                 :       1690 :   u = itor(muliu(R.Q,C/12), prec2);
     177                 :       1690 :   return rtor(mulrr(divri(u, R.T), sqrtr_abs(utor(C,prec2))), prec);
     178                 :            : }
     179                 :            : 
     180                 :            : #if 0 /* Much slower than binary splitting at least up to prec = 10^8 */
     181                 :            : /* Gauss - Brent-Salamin AGM iteration */
     182                 :            : static GEN
     183                 :            : pi_brent_salamin(long prec)
     184                 :            : {
     185                 :            :   GEN A, B, C;
     186                 :            :   pari_sp av2;
     187                 :            :   long i, G;
     188                 :            : 
     189                 :            :   G = - prec2nbits(prec);
     190                 :            :   incrprec(prec);
     191                 :            : 
     192                 :            :   A = real2n(-1, prec);
     193                 :            :   B = sqrtr_abs(A); /* = 1/sqrt(2) */
     194                 :            :   setexpo(A, 0);
     195                 :            :   C = real2n(-2, prec); av2 = avma;
     196                 :            :   for (i = 0;; i++)
     197                 :            :   {
     198                 :            :     GEN y, a, b, B_A = subrr(B, A);
     199                 :            :     pari_sp av3 = avma;
     200                 :            :     if (expo(B_A) < G) break;
     201                 :            :     a = addrr(A,B); shiftr_inplace(a, -1);
     202                 :            :     b = mulrr(A,B);
     203                 :            :     affrr(a, A);
     204                 :            :     affrr(sqrtr_abs(b), B); avma = av3;
     205                 :            :     y = sqrr(B_A); shiftr_inplace(y, i - 2);
     206                 :            :     affrr(subrr(C, y), C); avma = av2;
     207                 :            :   }
     208                 :            :   shiftr_inplace(C, 2);
     209                 :            :   return divrr(sqrr(addrr(A,B)), C);
     210                 :            : }
     211                 :            : GEN
     212                 :            : constpi(long prec)
     213                 :            : {
     214                 :            :   pari_sp av;
     215                 :            :   GEN tmp;
     216                 :            :   if (gpi && realprec(gpi) >= prec) return gpi;
     217                 :            : 
     218                 :            :   tmp = cgetr_block(prec);
     219                 :            :   av = avma;
     220                 :            :   affrr(pi_brent_salamin(prec), tmp);
     221                 :            :   swap_clone(&gpi, tmp);
     222                 :            :   avma = av;  return gpi;
     223                 :            : }
     224                 :            : #endif
     225                 :            : 
     226                 :            : GEN
     227                 :    5610889 : constpi(long prec)
     228                 :            : {
     229                 :            :   pari_sp av;
     230                 :            :   GEN tmp;
     231 [ +  + ][ +  + ]:    5610889 :   if (gpi && realprec(gpi) >= prec) return gpi;
     232                 :            : 
     233                 :       1690 :   av = avma;
     234                 :       1690 :   tmp = gclone(pi_ramanujan(prec));
     235                 :       1690 :   swap_clone(&gpi,tmp);
     236                 :    5610889 :   avma = av; return gpi;
     237                 :            : }
     238                 :            : 
     239                 :            : GEN
     240                 :    5610210 : mppi(long prec) { return rtor(constpi(prec), prec); }
     241                 :            : 
     242                 :            : /* Pi * 2^n */
     243                 :            : GEN
     244                 :    3804458 : Pi2n(long n, long prec)
     245                 :            : {
     246                 :    3804458 :   GEN x = mppi(prec); shiftr_inplace(x, n);
     247                 :    3804458 :   return x;
     248                 :            : }
     249                 :            : 
     250                 :            : /* I * Pi * 2^n */
     251                 :            : GEN
     252                 :       6886 : PiI2n(long n, long prec) { retmkcomplex(gen_0, Pi2n(n, prec)); }
     253                 :            : 
     254                 :            : /* 2I * Pi */
     255                 :            : GEN
     256                 :       4354 : PiI2(long prec) { return PiI2n(1, prec); }
     257                 :            : 
     258                 :            : /********************************************************************/
     259                 :            : /**                                                                **/
     260                 :            : /**                       EULER CONSTANT                           **/
     261                 :            : /**                                                                **/
     262                 :            : /********************************************************************/
     263                 :            : 
     264                 :            : GEN
     265                 :       7056 : consteuler(long prec)
     266                 :            : {
     267                 :            :   GEN u,v,a,b,tmpeuler;
     268                 :            :   long l, n1, n, k, x;
     269                 :            :   pari_sp av1, av2;
     270                 :            : 
     271 [ +  + ][ +  + ]:       7056 :   if (geuler && realprec(geuler) >= prec) return geuler;
     272                 :            : 
     273                 :        245 :   av1 = avma; tmpeuler = cgetr_block(prec);
     274                 :            : 
     275                 :        245 :   incrprec(prec);
     276                 :            : 
     277                 :        245 :   l = prec+EXTRAPRECWORD; x = (long) (1 + prec2nbits_mul(l, LOG2/4));
     278                 :        245 :   a = utor(x,l); u=logr_abs(a); setsigne(u,-1); affrr(u,a);
     279                 :        245 :   b = real_1(l);
     280                 :        245 :   v = real_1(l);
     281                 :        245 :   n = (long)(1+3.591*x); /* z=3.591: z*[ ln(z)-1 ]=1 */
     282                 :        245 :   n1 = minss(n, SQRTVERYBIGINT);
     283         [ +  - ]:        245 :   if (x < SQRTVERYBIGINT)
     284                 :            :   {
     285                 :        245 :     ulong xx = x*x;
     286                 :        245 :     av2 = avma;
     287         [ +  + ]:      77772 :     for (k=1; k<n1; k++)
     288                 :            :     {
     289                 :      77527 :       affrr(divru(mulur(xx,b),k*k), b);
     290                 :      77527 :       affrr(divru(addrr(divru(mulur(xx,a),k),b),k), a);
     291                 :      77527 :       affrr(addrr(u,a), u);
     292                 :      77527 :       affrr(addrr(v,b), v); avma = av2;
     293                 :            :     }
     294         [ +  + ]:        490 :     for (   ; k<=n; k++)
     295                 :            :     {
     296                 :        245 :       affrr(divru(divru(mulur(xx,b),k),k), b);
     297                 :        245 :       affrr(divru(addrr(divru(mulur(xx,a),k),b),k), a);
     298                 :        245 :       affrr(addrr(u,a), u);
     299                 :        245 :       affrr(addrr(v,b), v); avma = av2;
     300                 :            :     }
     301                 :            :   }
     302                 :            :   else
     303                 :            :   {
     304                 :          0 :     GEN xx = sqru(x);
     305                 :          0 :     av2 = avma;
     306         [ #  # ]:          0 :     for (k=1; k<n1; k++)
     307                 :            :     {
     308                 :          0 :       affrr(divru(mulir(xx,b),k*k), b);
     309                 :          0 :       affrr(divru(addrr(divru(mulir(xx,a),k),b),k), a);
     310                 :          0 :       affrr(addrr(u,a), u);
     311                 :          0 :       affrr(addrr(v,b), v); avma = av2;
     312                 :            :     }
     313         [ #  # ]:          0 :     for (   ; k<=n; k++)
     314                 :            :     {
     315                 :          0 :       affrr(divru(divru(mulir(xx,b),k),k), b);
     316                 :          0 :       affrr(divru(addrr(divru(mulir(xx,a),k),b),k), a);
     317                 :          0 :       affrr(addrr(u,a), u);
     318                 :          0 :       affrr(addrr(v,b), v); avma = av2;
     319                 :            :     }
     320                 :            :   }
     321                 :        245 :   divrrz(u,v,tmpeuler);
     322                 :        245 :   swap_clone(&geuler,tmpeuler);
     323                 :       7056 :   avma = av1; return geuler;
     324                 :            : }
     325                 :            : 
     326                 :            : GEN
     327                 :       7056 : mpeuler(long prec) { return rtor(consteuler(prec), prec); }
     328                 :            : 
     329                 :            : /********************************************************************/
     330                 :            : /**                                                                **/
     331                 :            : /**                       CATALAN CONSTANT                         **/
     332                 :            : /**                                                                **/
     333                 :            : /********************************************************************/
     334                 :            : /* 8G = 3\sum_{n>=0} 1/(binomial(2n,n)(2n+1)^2) + Pi log(2+sqrt(3)) */
     335                 :            : static GEN
     336                 :          7 : catalan(long prec)
     337                 :            : {
     338                 :          7 :   long i, nmax = prec2nbits(prec) >> 1;
     339                 :            :   struct abpq_res R;
     340                 :            :   struct abpq A;
     341                 :            :   GEN u, v;
     342                 :          7 :   abpq_init(&A, nmax);
     343                 :          7 :   A.a[0] = A.b[0] = A.p[0] = A.q[0] = gen_1;
     344         [ +  + ]:       5831 :   for (i = 1; i <= nmax; i++)
     345                 :            :   {
     346                 :       5824 :     A.a[i] = gen_1;
     347                 :       5824 :     A.b[i] = utoipos((i<<1)+1);
     348                 :       5824 :     A.p[i] = utoipos(i);
     349                 :       5824 :     A.q[i] = utoipos((i<<2)+2);
     350                 :            :   }
     351                 :          7 :   abpq_sum(&R, 0, nmax, &A);
     352                 :          7 :   u = mulur(3, rdivii(R.T, mulii(R.B,R.Q),prec));
     353                 :          7 :   v = mulrr(mppi(prec), logr_abs(addrs(sqrtr_abs(utor(3,prec)), 2)));
     354                 :          7 :   u = addrr(u, v); shiftr_inplace(u, -3);
     355                 :          7 :   return u;
     356                 :            : }
     357                 :            : 
     358                 :            : GEN
     359                 :          7 : constcatalan(long prec)
     360                 :            : {
     361                 :          7 :   pari_sp av = avma;
     362                 :            :   GEN tmp;
     363 [ -  + ][ #  # ]:          7 :   if (gcatalan && realprec(gcatalan) >= prec) return gcatalan;
     364                 :          7 :   tmp = gclone(catalan(prec));
     365                 :          7 :   swap_clone(&gcatalan,tmp);
     366                 :          7 :   avma = av; return gcatalan;
     367                 :            : }
     368                 :            : 
     369                 :            : GEN
     370                 :          7 : mpcatalan(long prec) { return rtor(constcatalan(prec), prec); }
     371                 :            : 
     372                 :            : /********************************************************************/
     373                 :            : /**                                                                **/
     374                 :            : /**          TYPE CONVERSION FOR TRANSCENDENTAL FUNCTIONS          **/
     375                 :            : /**                                                                **/
     376                 :            : /********************************************************************/
     377                 :            : static GEN
     378                 :      33627 : transvec(GEN (*f)(GEN,long), GEN x, long prec)
     379                 :            : {
     380                 :            :   long i, l;
     381                 :      33627 :   GEN y = cgetg_copy(x, &l);
     382         [ +  + ]:     141392 :   for (i=1; i<l; i++) gel(y,i) = f(gel(x,i),prec);
     383                 :      33627 :   return y;
     384                 :            : }
     385                 :            : 
     386                 :            : GEN
     387                 :     578543 : trans_eval(const char *fun, GEN (*f)(GEN,long), GEN x, long prec)
     388                 :            : {
     389                 :     578543 :   pari_sp av = avma;
     390         [ -  + ]:     578543 :   if (prec < 3) pari_err_BUG("trans_eval [prec < 3]");
     391   [ +  +  +  +  :     578543 :   switch(typ(x))
                   +  + ]
     392                 :            :   {
     393                 :     541955 :     case t_INT:    x = f(itor(x,prec),prec); break;
     394                 :       2905 :     case t_FRAC:   x = f(fractor(x, prec),prec); break;
     395                 :          7 :     case t_QUAD:   x = f(quadtofp(x,prec),prec); break;
     396                 :          7 :     case t_POLMOD: x = transvec(f, polmod_to_embed(x,prec), prec); break;
     397                 :            :     case t_VEC:
     398                 :            :     case t_COL:
     399                 :      33620 :     case t_MAT: return transvec(f, x, prec);
     400                 :         49 :     default: pari_err_TYPE(fun,x); return NULL;
     401                 :            :   }
     402                 :     578473 :   return gerepileupto(av, x);
     403                 :            : }
     404                 :            : 
     405                 :            : /*******************************************************************/
     406                 :            : /*                                                                 */
     407                 :            : /*                            POWERING                             */
     408                 :            : /*                                                                 */
     409                 :            : /*******************************************************************/
     410                 :            : /* x a t_REAL 0, return exp(x) */
     411                 :            : static GEN
     412                 :      31085 : mpexp0(GEN x)
     413                 :            : {
     414                 :      31085 :   long e = expo(x);
     415         [ -  + ]:      31085 :   return e >= 0? real_0_bit(e): real_1(nbits2prec(-e));
     416                 :            : }
     417                 :            : static GEN
     418                 :       1729 : powr0(GEN x)
     419         [ +  - ]:       1729 : { return signe(x)? real_1(realprec(x)): mpexp0(x); }
     420                 :            : 
     421                 :            : /* x t_POL or t_SER, return scalarpol(RgX_get_1(x)) */
     422                 :            : static GEN
     423                 :     106057 : scalarpol_get_1(GEN x)
     424                 :            : {
     425                 :     106057 :   GEN y = cgetg(3,t_POL);
     426                 :     106057 :   y[1] = evalvarn(varn(x)) | evalsigne(1);
     427                 :     106057 :   gel(y,2) = RgX_get_1(x); return y;
     428                 :            : }
     429                 :            : /* to be called by the generic function gpowgs(x,s) when s = 0 */
     430                 :            : static GEN
     431                 :     127288 : gpowg0(GEN x)
     432                 :            : {
     433                 :            :   long lx, i;
     434                 :            :   GEN y;
     435                 :            : 
     436   [ +  +  +  +  :     127288 :   switch(typ(x))
          +  +  +  +  +  
             +  +  +  + ]
     437                 :            :   {
     438                 :            :     case t_INT: case t_REAL: case t_FRAC: case t_PADIC:
     439                 :      13279 :       return gen_1;
     440                 :            : 
     441                 :          7 :     case t_QUAD: x++; /*fall through*/
     442                 :            :     case t_COMPLEX: {
     443                 :       3206 :       pari_sp av = avma;
     444                 :       3206 :       GEN a = gpowg0(gel(x,1));
     445                 :       3206 :       GEN b = gpowg0(gel(x,2));
     446         [ +  + ]:       3206 :       if (a == gen_1) return b;
     447         [ +  + ]:         14 :       if (b == gen_1) return a;
     448                 :          7 :       return gerepileupto(av, gmul(a,b));
     449                 :            :     }
     450                 :            :     case t_INTMOD:
     451                 :         56 :       y = cgetg(3,t_INTMOD);
     452                 :         56 :       gel(y,1) = icopy(gel(x,1));
     453                 :         56 :       gel(y,2) = gen_1; return y;
     454                 :            : 
     455                 :       4396 :     case t_FFELT: return FF_1(x);
     456                 :            : 
     457                 :            :     case t_POLMOD:
     458                 :         28 :       y = cgetg(3,t_POLMOD);
     459                 :         28 :       gel(y,1) = gcopy(gel(x,1));
     460                 :         28 :       gel(y,2) = scalarpol_get_1(gel(x,1)); return y;
     461                 :            : 
     462                 :            :     case t_RFRAC:
     463                 :          7 :       return scalarpol_get_1(gel(x,2));
     464                 :            :     case t_POL: case t_SER:
     465                 :     106022 :       return scalarpol_get_1(x);
     466                 :            : 
     467                 :            :     case t_MAT:
     468         [ +  + ]:         14 :       lx=lg(x); if (lx==1) return cgetg(1,t_MAT);
     469         [ -  + ]:          7 :       if (lx != lgcols(x)) pari_err_DIM("gpow");
     470                 :          7 :       y = matid(lx-1);
     471         [ +  + ]:         14 :       for (i=1; i<lx; i++) gcoeff(y,i,i) = gpowg0(gcoeff(x,i,i));
     472                 :          7 :       return y;
     473                 :         14 :     case t_QFR: return qfr_1(x);
     474                 :        252 :     case t_QFI: return qfi_1(x);
     475                 :          7 :     case t_VECSMALL: return identity_perm(lg(x) - 1);
     476                 :            :   }
     477                 :          7 :   pari_err_TYPE("gpow",x);
     478                 :     127281 :   return NULL; /* not reached */
     479                 :            : }
     480                 :            : 
     481                 :            : static GEN
     482                 :    3878554 : _sqr(void *data /* ignored */, GEN x) { (void)data; return gsqr(x); }
     483                 :            : static GEN
     484                 :    1473180 : _mul(void *data /* ignored */, GEN x, GEN y) { (void)data; return gmul(x,y); }
     485                 :            : static GEN
     486                 :        280 : _one(void *x) { return gpowg0((GEN) x); }
     487                 :            : static GEN
     488                 :    3948047 : _sqri(void *data /* ignored */, GEN x) { (void)data; return sqri(x); }
     489                 :            : static GEN
     490                 :    2603189 : _muli(void *data /* ignored */, GEN x, GEN y) { (void)data; return mulii(x,y); }
     491                 :            : static GEN
     492                 :   26359626 : _sqrr(void *data /* ignored */, GEN x) { (void)data; return sqrr(x); }
     493                 :            : static GEN
     494                 :     562528 : _mulr(void *data /* ignored */, GEN x, GEN y) { (void)data; return mulrr(x,y); }
     495                 :            : static GEN
     496                 :      16255 : _oner(void *data /* prec */) { return real_1( *(long*) data); }
     497                 :            : 
     498                 :            : /* INTEGER POWERING (a^n for integer a != 0 and integer n > 0)
     499                 :            :  *
     500                 :            :  * Use left shift binary algorithm (RS is wasteful: multiplies big numbers,
     501                 :            :  * with LS one of them is the base, hence small). Sign of result is set
     502                 :            :  * to s (= 1,-1). Makes life easier for caller, which otherwise might do a
     503                 :            :  * setsigne(gen_1 / gen_m1) */
     504                 :            : static GEN
     505                 :   15746610 : powiu_sign(GEN a, ulong N, long s)
     506                 :            : {
     507                 :            :   pari_sp av;
     508                 :            :   GEN y;
     509                 :            : 
     510         [ +  + ]:   15746610 :   if (lgefint(a) == 3)
     511                 :            :   { /* easy if |a| < 3 */
     512                 :   12779434 :     ulong q = a[2];
     513 [ +  + ][ +  + ]:   12779434 :     if (q == 1) return (s>0)? gen_1: gen_m1;
     514         [ +  + ]:    4771849 :     if (q == 2) { a = int2u(N); setsigne(a,s); return a; }
     515                 :    3448860 :     q = upowuu(q, N);
     516 [ +  + ][ +  + ]:    3448860 :     if (q) return s>0? utoipos(q): utoineg(q);
     517                 :            :   }
     518         [ +  + ]:    3480544 :   if (N <= 2) {
     519         [ +  + ]:    2243852 :     if (N == 2) return sqri(a);
     520                 :       3382 :     a = icopy(a); setsigne(a,s); return a;
     521                 :            :   }
     522                 :    1236692 :   av = avma;
     523                 :    1236692 :   y = gen_powu_i(a, N, NULL, &_sqri, &_muli);
     524                 :   15746610 :   setsigne(y,s); return gerepileuptoint(av, y);
     525                 :            : }
     526                 :            : /* a^n */
     527                 :            : GEN
     528                 :   15279452 : powiu(GEN a, ulong n)
     529                 :            : {
     530                 :            :   long s;
     531         [ +  + ]:   15279452 :   if (!n) return gen_1;
     532                 :   15000379 :   s = signe(a);
     533         [ +  + ]:   15000379 :   if (!s) return gen_0;
     534 [ +  + ][ +  + ]:   15279452 :   return powiu_sign(a, n, (s < 0 && odd(n))? -1: 1);
     535                 :            : }
     536                 :            : GEN
     537                 :   11665936 : powis(GEN a, long n)
     538                 :            : {
     539                 :            :   long s;
     540                 :            :   GEN t, y;
     541         [ +  + ]:   11665936 :   if (n >= 0) return powiu(a, n);
     542                 :      68609 :   s = signe(a);
     543         [ -  + ]:      68609 :   if (!s) pari_err_INV("powis",gen_0);
     544 [ +  + ][ -  + ]:      68609 :   t = (s < 0 && odd(n))? gen_m1: gen_1;
     545         [ +  + ]:      68609 :   if (is_pm1(a)) return t;
     546                 :            :   /* n < 0, |a| > 1 */
     547                 :      68497 :   y = cgetg(3,t_FRAC);
     548                 :      68497 :   gel(y,1) = t;
     549                 :      68497 :   gel(y,2) = powiu_sign(a, -n, 1); /* force denominator > 0 */
     550                 :   11665936 :   return y;
     551                 :            : }
     552                 :            : GEN
     553                 :    4052922 : powuu(ulong p, ulong N)
     554                 :            : {
     555                 :    4052922 :   pari_sp av = avma;
     556                 :    4052922 :   long P[] = {evaltyp(t_INT)|_evallg(3), evalsigne(1)|evallgefint(3),0};
     557                 :            :   ulong pN;
     558                 :            :   GEN y;
     559         [ +  + ]:    4052922 :   if (N <= 2)
     560                 :            :   {
     561         [ +  + ]:    2095994 :     if (N == 2) return sqru(p);
     562         [ +  + ]:    1587779 :     if (N == 1) return utoipos(p);
     563                 :     815325 :     return gen_1;
     564                 :            :   }
     565         [ -  + ]:    1956928 :   if (!p) return gen_0;
     566                 :    1956928 :   pN = upowuu(p, N);
     567         [ +  + ]:    1956928 :   if (pN) return utoipos(pN);
     568         [ +  + ]:     645108 :   if (p == 2) return int2u(N);
     569                 :     644409 :   P[2] = p; av = avma;
     570                 :     644409 :   y = gen_powu_i(P, N, NULL, &_sqri, &_muli);
     571                 :    4052922 :   return gerepileuptoint(av, y);
     572                 :            : }
     573                 :            : 
     574                 :            : /* return 0 if overflow */
     575                 :            : static ulong
     576         [ +  + ]:    1658972 : usqru(ulong p) { return p & HIGHMASK? 0: p*p; }
     577                 :            : ulong
     578                 :    7820951 : upowuu(ulong p, ulong k)
     579                 :            : {
     580                 :            : #ifdef LONG_IS_64BIT
     581                 :    6768852 :   const ulong CUTOFF3 = 2642245;
     582                 :    6768852 :   const ulong CUTOFF4 = 65535;
     583                 :    6768852 :   const ulong CUTOFF5 = 7131;
     584                 :    6768852 :   const ulong CUTOFF6 = 1625;
     585                 :    6768852 :   const ulong CUTOFF7 = 565;
     586                 :    6768852 :   const ulong CUTOFF8 = 255;
     587                 :    6768852 :   const ulong CUTOFF9 = 138;
     588                 :    6768852 :   const ulong CUTOFF10 = 84;
     589                 :    6768852 :   const ulong CUTOFF11 = 56;
     590                 :    6768852 :   const ulong CUTOFF12 = 40;
     591                 :    6768852 :   const ulong CUTOFF13 = 30;
     592                 :    6768852 :   const ulong CUTOFF14 = 23;
     593                 :    6768852 :   const ulong CUTOFF15 = 19;
     594                 :    6768852 :   const ulong CUTOFF16 = 15;
     595                 :    6768852 :   const ulong CUTOFF17 = 13;
     596                 :    6768852 :   const ulong CUTOFF18 = 11;
     597                 :    6768852 :   const ulong CUTOFF19 = 10;
     598                 :    6768852 :   const ulong CUTOFF20 =  9;
     599                 :            : #else
     600                 :    1052099 :   const ulong CUTOFF3 = 1625;
     601                 :    1052099 :   const ulong CUTOFF4 =  255;
     602                 :    1052099 :   const ulong CUTOFF5 =   84;
     603                 :    1052099 :   const ulong CUTOFF6 =   40;
     604                 :    1052099 :   const ulong CUTOFF7 =   23;
     605                 :    1052099 :   const ulong CUTOFF8 =   15;
     606                 :    1052099 :   const ulong CUTOFF9 =   11;
     607                 :    1052099 :   const ulong CUTOFF10 =   9;
     608                 :    1052099 :   const ulong CUTOFF11 =   7;
     609                 :    1052099 :   const ulong CUTOFF12 =   6;
     610                 :    1052099 :   const ulong CUTOFF13 =   5;
     611                 :    1052099 :   const ulong CUTOFF14 =   4;
     612                 :    1052099 :   const ulong CUTOFF15 =   4;
     613                 :    1052099 :   const ulong CUTOFF16 =   3;
     614                 :    1052099 :   const ulong CUTOFF17 =   3;
     615                 :    1052099 :   const ulong CUTOFF18 =   3;
     616                 :    1052099 :   const ulong CUTOFF19 =   3;
     617                 :    1052099 :   const ulong CUTOFF20 =   3;
     618                 :            : #endif
     619                 :            : 
     620         [ +  + ]:    7820951 :   if (p <= 2)
     621                 :            :   {
     622         [ +  + ]:     696963 :     if (p < 2) return p;
     623         [ +  + ]:     694282 :     return k < BITS_IN_LONG? 1UL<<k: 0;
     624                 :            :   }
     625   [ +  +  +  +  :    7123988 :   switch(k)
          +  +  +  +  +  
          +  +  +  +  +  
          +  +  +  +  +  
                +  +  + ]
     626                 :            :   {
     627                 :            :     ulong p2, p3, p4, p5, p8;
     628                 :     554333 :     case 0:  return 1;
     629                 :    1517292 :     case 1:  return p;
     630                 :    1658972 :     case 2:  return usqru(p);
     631         [ +  + ]:     908119 :     case 3:  if (p > CUTOFF3) return 0; return p*p*p;
     632         [ +  + ]:     652953 :     case 4:  if (p > CUTOFF4) return 0; p2=p*p; return p2*p2;
     633         [ +  + ]:     280717 :     case 5:  if (p > CUTOFF5) return 0; p2=p*p; return p2*p2*p;
     634         [ +  + ]:     496485 :     case 6:  if (p > CUTOFF6) return 0; p2=p*p; return p2*p2*p2;
     635         [ +  + ]:      81265 :     case 7:  if (p > CUTOFF7) return 0; p2=p*p; return p2*p2*p2*p;
     636         [ +  + ]:      86618 :     case 8:  if (p > CUTOFF8) return 0; p2=p*p; p4=p2*p2; return p4*p4;
     637         [ +  + ]:      51911 :     case 9:  if (p > CUTOFF9) return 0; p2=p*p; p4=p2*p2; return p4*p4*p;
     638         [ +  + ]:      54765 :     case 10: if (p > CUTOFF10)return 0; p2=p*p; p4=p2*p2; return p4*p4*p2;
     639         [ +  + ]:      23215 :     case 11: if (p > CUTOFF11)return 0; p2=p*p; p4=p2*p2; return p4*p4*p2*p;
     640         [ +  + ]:      30562 :     case 12: if (p > CUTOFF12)return 0; p2=p*p; p4=p2*p2; return p4*p4*p4;
     641         [ +  + ]:      67802 :     case 13: if (p > CUTOFF13)return 0; p2=p*p; p4=p2*p2; return p4*p4*p4*p;
     642         [ +  + ]:      58317 :     case 14: if (p > CUTOFF14)return 0; p2=p*p; p4=p2*p2; return p4*p4*p4*p2;
     643         [ +  + ]:      66699 :     case 15: if (p > CUTOFF15)return 0;
     644                 :       9619 :       p2=p*p; p3=p2*p; p5=p3*p2; return p5*p5*p5;
     645         [ +  + ]:      53404 :     case 16: if (p > CUTOFF16)return 0;
     646                 :       9599 :       p2=p*p; p4=p2*p2; p8=p4*p4; return p8*p8;
     647         [ +  + ]:      32048 :     case 17: if (p > CUTOFF17)return 0;
     648                 :       9665 :       p2=p*p; p4=p2*p2; p8=p4*p4; return p*p8*p8;
     649         [ +  + ]:      22238 :     case 18: if (p > CUTOFF18)return 0;
     650                 :       6665 :       p2=p*p; p4=p2*p2; p8=p4*p4; return p2*p8*p8;
     651         [ +  + ]:      18820 :     case 19: if (p > CUTOFF19)return 0;
     652                 :       7388 :       p2=p*p; p4=p2*p2; p8=p4*p4; return p*p2*p8*p8;
     653         [ +  + ]:      11598 :     case 20: if (p > CUTOFF20)return 0;
     654                 :       3420 :       p2=p*p; p4=p2*p2; p8=p4*p4; return p4*p8*p8;
     655                 :            :   }
     656                 :            : #ifdef LONG_IS_64BIT
     657   [ +  +  +  +  :     337629 :   switch(p)
                   +  + ]
     658                 :            :   {
     659         [ +  + ]:      23388 :     case 3: if (k > 40) return 0;
     660                 :      16656 :       break;
     661         [ +  + ]:        756 :     case 4: if (k > 31) return 0;
     662                 :        150 :       return 1UL<<(2*k);
     663         [ +  + ]:       6042 :     case 5: if (k > 27) return 0;
     664                 :       2526 :       break;
     665         [ +  + ]:        642 :     case 6: if (k > 24) return 0;
     666                 :         24 :       break;
     667         [ +  + ]:       9414 :     case 7: if (k > 22) return 0;
     668                 :       1020 :       break;
     669                 :     297387 :     default: return 0;
     670                 :            :   }
     671                 :            :   /* no overflow */
     672                 :            :   {
     673                 :      20226 :     ulong q = upowuu(p, k >> 1);
     674                 :      20226 :     q *= q ;
     675         [ +  + ]:    6768852 :     return odd(k)? q*p: q;
     676                 :            :   }
     677                 :            : #else
     678                 :    1052099 :   return 0;
     679                 :            : #endif
     680                 :            : }
     681                 :            : 
     682                 :            : typedef struct {
     683                 :            :   long prec, a;
     684                 :            :   GEN (*sqr)(GEN);
     685                 :            :   GEN (*mulug)(ulong,GEN);
     686                 :            : } sr_muldata;
     687                 :            : 
     688                 :            : static GEN
     689                 :     322615 : _rpowuu_sqr(void *data, GEN x)
     690                 :            : {
     691                 :     322615 :   sr_muldata *D = (sr_muldata *)data;
     692 [ +  + ][ +  + ]:     322615 :   if (typ(x) == t_INT && lgefint(x) >= D->prec)
     693                 :            :   { /* switch to t_REAL */
     694                 :       3718 :     D->sqr   = &sqrr;
     695                 :       3718 :     D->mulug = &mulur; x = itor(x, D->prec);
     696                 :            :   }
     697                 :     322615 :   return D->sqr(x);
     698                 :            : }
     699                 :            : 
     700                 :            : static GEN
     701                 :     103618 : _rpowuu_msqr(void *data, GEN x)
     702                 :            : {
     703                 :     103618 :   GEN x2 = _rpowuu_sqr(data, x);
     704                 :     103618 :   sr_muldata *D = (sr_muldata *)data;
     705                 :     103618 :   return D->mulug(D->a, x2);
     706                 :            : }
     707                 :            : 
     708                 :            : /* return a^n as a t_REAL of precision prec. Assume a > 0, n > 0 */
     709                 :            : GEN
     710                 :     110060 : rpowuu(ulong a, ulong n, long prec)
     711                 :            : {
     712                 :            :   pari_sp av;
     713                 :            :   GEN y, z;
     714                 :            :   sr_muldata D;
     715                 :            : 
     716         [ -  + ]:     110060 :   if (a == 1) return real_1(prec);
     717         [ -  + ]:     110060 :   if (a == 2) return real2n(n, prec);
     718         [ +  + ]:     110060 :   if (n == 1) return utor(a, prec);
     719                 :     109402 :   z = cgetr(prec);
     720                 :     109402 :   av = avma;
     721                 :     109402 :   D.sqr   = &sqri;
     722                 :     109402 :   D.mulug = &mului;
     723                 :     109402 :   D.prec = prec;
     724                 :     109402 :   D.a = (long)a;
     725                 :     109402 :   y = gen_powu_fold_i(utoipos(a), n, (void*)&D, &_rpowuu_sqr, &_rpowuu_msqr);
     726                 :     110060 :   mpaff(y, z); avma = av; return z;
     727                 :            : }
     728                 :            : 
     729                 :            : GEN
     730                 :   25578702 : powrs(GEN x, long n)
     731                 :            : {
     732                 :   25578702 :   pari_sp av = avma;
     733                 :            :   GEN y;
     734         [ -  + ]:   25578702 :   if (!n) return powr0(x);
     735                 :   25578702 :   y = gen_powu_i(x, (ulong)labs(n), NULL, &_sqrr, &_mulr);
     736         [ +  + ]:   25578702 :   if (n < 0) y = invr(y);
     737                 :   25578702 :   return gerepileuptoleaf(av,y);
     738                 :            : }
     739                 :            : GEN
     740                 :     144710 : powru(GEN x, ulong n)
     741                 :            : {
     742                 :     144710 :   pari_sp av = avma;
     743                 :            :   GEN y;
     744         [ +  + ]:     144710 :   if (!n) return powr0(x);
     745                 :     143499 :   y = gen_powu_i(x, n, NULL, &_sqrr, &_mulr);
     746                 :     144710 :   return gerepileuptoleaf(av,y);
     747                 :            : }
     748                 :            : 
     749                 :            : GEN
     750                 :      16255 : powersr(GEN x, long n)
     751                 :            : {
     752                 :      16255 :   long prec = realprec(x);
     753                 :      16255 :   return gen_powers(x, n, 1, &prec, &_sqrr, &_mulr, &_oner);
     754                 :            : }
     755                 :            : 
     756                 :            : /* x^(s/2), assume x t_REAL */
     757                 :            : GEN
     758                 :          0 : powrshalf(GEN x, long s)
     759                 :            : {
     760         [ #  # ]:          0 :   if (s & 1) return sqrtr(powrs(x, s));
     761                 :          0 :   return powrs(x, s>>1);
     762                 :            : }
     763                 :            : /* x^(s/2), assume x t_REAL */
     764                 :            : GEN
     765                 :      11886 : powruhalf(GEN x, ulong s)
     766                 :            : {
     767         [ +  + ]:      11886 :   if (s & 1) return sqrtr(powru(x, s));
     768                 :      11886 :   return powru(x, s>>1);
     769                 :            : }
     770                 :            : /* x^(n/d), assume x t_REAL, return t_REAL */
     771                 :            : GEN
     772                 :        539 : powrfrac(GEN x, long n, long d)
     773                 :            : {
     774                 :            :   long z;
     775         [ +  + ]:        539 :   if (!n) return powr0(x);
     776         [ +  + ]:         21 :   z = cgcd(n, d); if (z > 1) { n /= z; d /= z; }
     777         [ +  + ]:         21 :   if (d == 1) return powrs(x, n);
     778                 :         14 :   x = powrs(x, n);
     779         [ -  + ]:         14 :   if (d == 2) return sqrtr(x);
     780                 :        539 :   return sqrtnr(x, d);
     781                 :            : }
     782                 :            : 
     783                 :            : /* assume x != 0 */
     784                 :            : static GEN
     785                 :     104132 : pow_monome(GEN x, long n)
     786                 :            : {
     787                 :     104132 :   long i, d, dx = degpol(x);
     788                 :            :   GEN A, b, y;
     789                 :            : 
     790         [ +  + ]:     104132 :   if (n < 0) { n = -n; y = cgetg(3, t_RFRAC); } else y = NULL;
     791                 :            : 
     792 [ +  + ][ +  + ]:     104132 :   if (HIGHWORD(dx) || HIGHWORD(n))
     793                 :          8 :   {
     794                 :            :     LOCAL_HIREMAINDER;
     795                 :          9 :     d = (long)mulll((ulong)dx, (ulong)n);
     796 [ +  + ][ -  + ]:          9 :     if (hiremainder || (d &~ LGBITS)) d = LGBITS; /* overflow */
         [ -  + ][ #  # ]
     797                 :          9 :     d += 2;
     798                 :            :   }
     799                 :            :   else
     800                 :     104123 :     d = dx*n + 2;
     801         [ +  + ]:     104132 :   if ((d + 1) & ~LGBITS) pari_err(e_OVERFLOW,"pow_monome [degree]");
     802                 :     104125 :   A = cgetg(d+1, t_POL); A[1] = x[1];
     803         [ +  + ]:    2681294 :   for (i=2; i < d; i++) gel(A,i) = gen_0;
     804                 :     104125 :   b = gpowgs(gel(x,dx+2), n); /* not memory clean if (n < 0) */
     805         [ +  + ]:     104125 :   if (!y) y = A;
     806                 :            :   else {
     807                 :      20433 :     GEN c = denom(b);
     808         [ -  + ]:      20433 :     gel(y,1) = c; if (c != gen_1) b = gmul(b,c);
     809                 :      20433 :     gel(y,2) = A;
     810                 :            :   }
     811                 :     104125 :   gel(A,d) = b; return y;
     812                 :            : }
     813                 :            : 
     814                 :            : /* x t_PADIC */
     815                 :            : static GEN
     816                 :       1358 : powps(GEN x, long n)
     817                 :            : {
     818                 :       1358 :   long e = n*valp(x), v;
     819                 :       1358 :   GEN t, y, mod, p = gel(x,2);
     820                 :            :   pari_sp av;
     821                 :            : 
     822         [ +  + ]:       1358 :   if (!signe(gel(x,4))) {
     823         [ +  + ]:         84 :     if (n < 0) pari_err_INV("powps",x);
     824                 :         77 :     return zeropadic(p, e);
     825                 :            :   }
     826                 :       1274 :   v = z_pval(n, p);
     827                 :            : 
     828                 :       1274 :   y = cgetg(5,t_PADIC);
     829                 :       1274 :   mod = gel(x,3);
     830         [ +  + ]:       1274 :   if (v == 0) mod = icopy(mod);
     831                 :            :   else
     832                 :            :   {
     833 [ +  + ][ +  + ]:        476 :     if (precp(x) == 1 && equaliu(p, 2)) v++;
     834                 :        476 :     mod = mulii(mod, powiu(p,v));
     835                 :        476 :     mod = gerepileuptoint((pari_sp)y, mod);
     836                 :            :   }
     837                 :       1274 :   y[1] = evalprecp(precp(x) + v) | evalvalp(e);
     838                 :       1274 :   gel(y,2) = icopy(p);
     839                 :       1274 :   gel(y,3) = mod;
     840                 :            : 
     841                 :       1274 :   av = avma; t = gel(x,4);
     842         [ +  + ]:       1274 :   if (n < 0) { t = Fp_inv(t, mod); n = -n; }
     843                 :       1274 :   t = Fp_powu(t, n, mod);
     844                 :       1274 :   gel(y,4) = gerepileuptoint(av, t);
     845                 :       1351 :   return y;
     846                 :            : }
     847                 :            : /* x t_PADIC */
     848                 :            : static GEN
     849                 :        161 : powp(GEN x, GEN n)
     850                 :            : {
     851                 :            :   long v;
     852                 :        161 :   GEN y, mod, p = gel(x,2);
     853                 :            : 
     854         [ -  + ]:        161 :   if (valp(x)) pari_err_OVERFLOW("valp()");
     855                 :            : 
     856         [ +  + ]:        161 :   if (!signe(gel(x,4))) {
     857         [ +  + ]:         14 :     if (signe(n) < 0) pari_err_INV("powp",x);
     858                 :          7 :     return zeropadic(p, 0);
     859                 :            :   }
     860                 :        147 :   v = Z_pval(n, p);
     861                 :            : 
     862                 :        147 :   y = cgetg(5,t_PADIC);
     863                 :        147 :   mod = gel(x,3);
     864         [ +  + ]:        147 :   if (v == 0) mod = icopy(mod);
     865                 :            :   else
     866                 :            :   {
     867                 :         70 :     mod = mulii(mod, powiu(p,v));
     868                 :         70 :     mod = gerepileuptoint((pari_sp)y, mod);
     869                 :            :   }
     870                 :        147 :   y[1] = evalprecp(precp(x) + v) | _evalvalp(0);
     871                 :        147 :   gel(y,2) = icopy(p);
     872                 :        147 :   gel(y,3) = mod;
     873                 :        147 :   gel(y,4) = Fp_pow(gel(x,4), n, mod);
     874                 :        154 :   return y;
     875                 :            : }
     876                 :            : static GEN
     877                 :      10557 : pow_polmod(GEN x, GEN n)
     878                 :            : {
     879                 :      10557 :   GEN z = cgetg(3, t_POLMOD), a = gel(x,2), T = gel(x,1);
     880                 :      10557 :   gel(z,1) = gcopy(T);
     881 [ +  + ][ +  - ]:      10557 :   if (typ(a) != t_POL || varn(a) != varn(T) || lg(a) <= 3)
                 [ +  + ]
     882                 :       1120 :     a = powgi(a, n);
     883                 :            :   else {
     884                 :       9437 :     pari_sp av = avma;
     885                 :       9437 :     GEN p = NULL;
     886 [ +  + ][ +  + ]:       9437 :     if (RgX_is_FpX(T, &p) && RgX_is_FpX(a, &p) && p)
                 [ +  + ]
     887                 :            :     {
     888                 :       7602 :       T = RgX_to_FpX(T, p); a = RgX_to_FpX(a, p);
     889         [ +  + ]:       7602 :       if (lgefint(p) == 3)
     890                 :            :       {
     891                 :       7595 :         ulong pp = p[2];
     892                 :       7595 :         a = Flxq_pow(ZX_to_Flx(a, pp), n, ZX_to_Flx(T, pp), pp);
     893                 :       7595 :         a = Flx_to_ZX(a);
     894                 :            :       }
     895                 :            :       else
     896                 :          7 :         a = FpXQ_pow(a, n, T, p);
     897                 :       7602 :       a = FpX_to_mod(a, p);
     898                 :       7602 :       a = gerepileupto(av, a);
     899                 :            :     }
     900                 :            :     else
     901                 :            :     {
     902                 :       1835 :       avma = av;
     903                 :       9437 :       a = RgXQ_pow(a, n, gel(z,1));
     904                 :            :     }
     905                 :            :   }
     906                 :      10557 :   gel(z,2) = a; return z;
     907                 :            : }
     908                 :            : 
     909                 :            : GEN
     910                 :   40294030 : gpowgs(GEN x, long n)
     911                 :            : {
     912                 :            :   long m;
     913                 :            :   pari_sp av;
     914                 :            :   GEN y;
     915                 :            : 
     916         [ +  + ]:   40294030 :   if (n == 0) return gpowg0(x);
     917         [ +  + ]:   40173441 :   if (n == 1)
     918      [ +  +  + ]:    1384840 :     switch (typ(x)) {
     919                 :       2415 :       case t_QFI: return redimag(x);
     920                 :         14 :       case t_QFR: return redreal(x);
     921                 :    1382411 :       default: return gcopy(x);
     922                 :            :     }
     923         [ +  + ]:   38788601 :   if (n ==-1) return ginv(x);
     924   [ +  +  +  +  :   38738537 :   switch(typ(x))
             +  +  +  +  
                      + ]
     925                 :            :   {
     926                 :   11623556 :     case t_INT: return powis(x,n);
     927                 :   25578177 :     case t_REAL: return powrs(x,n);
     928                 :            :     case t_INTMOD:
     929                 :      14336 :       y = cgetg(3,t_INTMOD); gel(y,1) = icopy(gel(x,1));
     930                 :      14336 :       gel(y,2) = Fp_pows(gel(x,2), n, gel(x,1));
     931                 :      14336 :       return y;
     932                 :            :     case t_FRAC:
     933                 :            :     {
     934                 :     344435 :       GEN a = gel(x,1), b = gel(x,2);
     935 [ +  + ][ +  + ]:     344435 :       long s = (signe(a) < 0 && odd(n))? -1: 1;
     936         [ +  + ]:     344435 :       if (n < 0) {
     937                 :         14 :         n = -n;
     938         [ +  + ]:         14 :         if (is_pm1(a)) return powiu_sign(b, n, s); /* +-1/x[2] inverts to t_INT */
     939                 :          7 :         swap(a, b);
     940                 :            :       }
     941                 :     344428 :       y = cgetg(3, t_FRAC);
     942                 :     344428 :       gel(y,1) = powiu_sign(a, n, s);
     943                 :     344428 :       gel(y,2) = powiu_sign(b, n, 1);
     944                 :     344428 :       return y;
     945                 :            :     }
     946                 :       1358 :     case t_PADIC: return powps(x, n);
     947                 :            :     case t_RFRAC:
     948                 :            :     {
     949                 :         21 :       av = avma; y = cgetg(3, t_RFRAC); m = labs(n);
     950                 :         21 :       gel(y,1) = gpowgs(gel(x,1),m);
     951                 :         21 :       gel(y,2) = gpowgs(gel(x,2),m);
     952         [ -  + ]:         21 :       if (n < 0) y = ginv(y);
     953                 :         21 :       return gerepileupto(av,y);
     954                 :            :     }
     955                 :            :     case t_POLMOD: {
     956                 :      10550 :       long N[] = {evaltyp(t_INT) | _evallg(3),0,0};
     957                 :      10550 :       affsi(n,N); return pow_polmod(x, N);
     958                 :            :     }
     959                 :            :     case t_POL:
     960         [ +  + ]:     128520 :       if (RgX_is_monomial(x)) return pow_monome(x, n);
     961                 :            :     default: {
     962                 :    1061972 :       pari_sp av = avma;
     963                 :    1061972 :       y = gen_powu_i(x, (ulong)labs(n), NULL, &_sqr, &_mul);
     964         [ +  + ]:    1061972 :       if (n < 0) y = ginv(y);
     965                 :   40294002 :       return gerepileupto(av,y);
     966                 :            :     }
     967                 :            :   }
     968                 :            : }
     969                 :            : 
     970                 :            : /* n a t_INT */
     971                 :            : GEN
     972                 :    9666458 : powgi(GEN x, GEN n)
     973                 :            : {
     974                 :            :   GEN y;
     975                 :            : 
     976         [ +  + ]:    9666458 :   if (!is_bigint(n)) return gpowgs(x, itos(n));
     977                 :            :   /* probable overflow for non-modular types (typical exception: (X^0)^N) */
     978   [ +  +  +  +  :        285 :   switch(typ(x))
             +  +  +  + ]
     979                 :            :   {
     980                 :            :     case t_INTMOD:
     981                 :         28 :       y = cgetg(3,t_INTMOD); gel(y,1) = icopy(gel(x,1));
     982                 :         28 :       gel(y,2) = Fp_pow(gel(x,2), n, gel(x,1));
     983                 :         28 :       return y;
     984                 :         26 :     case t_FFELT: return FF_pow(x,n);
     985                 :        161 :     case t_PADIC: return powp(x, n);
     986                 :            : 
     987                 :            :     case t_INT:
     988 [ +  + ][ +  + ]:         35 :       if (is_pm1(x)) return (signe(x) < 0 && mpodd(n))? gen_m1: gen_1;
                 [ +  + ]
     989         [ +  + ]:         14 :       if (signe(x)) pari_err_OVERFLOW("lg()");
     990         [ -  + ]:          7 :       if (signe(n) < 0) pari_err_INV("powgi",gen_0);
     991                 :          7 :       return gen_0;
     992                 :            :     case t_FRAC:
     993                 :          7 :       pari_err_OVERFLOW("lg()");
     994                 :            : 
     995                 :         14 :     case t_QFR: return qfrpow(x, n);
     996                 :          7 :     case t_POLMOD: return pow_polmod(x, n);
     997                 :            :     default: {
     998                 :          7 :       pari_sp av = avma;
     999                 :          7 :       y = gen_pow(x, n, NULL, &_sqr, &_mul);
    1000         [ -  + ]:          7 :       if (signe(n) < 0) y = ginv(y);
    1001                 :    9666416 :       return gerepileupto(av,y);
    1002                 :            :     }
    1003                 :            :   }
    1004                 :            : }
    1005                 :            : 
    1006                 :            : /* Assume x = 1 + O(t), n a scalar. Return x^n */
    1007                 :            : static GEN
    1008                 :        161 : ser_pow_1(GEN x, GEN n)
    1009                 :            : {
    1010                 :            :   long lx, mi, i, j, d;
    1011                 :        161 :   GEN y = cgetg_copy(x, &lx), X = x+2, Y = y + 2;
    1012                 :        161 :   y[1] = evalsigne(1) | _evalvalp(0) | evalvarn(varn(x));
    1013 [ +  + ][ +  + ]:        189 :   d = mi = lx-3; while (mi>=1 && isrationalzero(gel(X,mi))) mi--;
    1014                 :        161 :   gel(Y,0) = gen_1;
    1015         [ +  + ]:       1883 :   for (i=1; i<=d; i++)
    1016                 :            :   {
    1017                 :       1722 :     pari_sp av = avma;
    1018                 :       1722 :     GEN s = gen_0;
    1019         [ +  + ]:      17192 :     for (j=1; j<=minss(i,mi); j++)
    1020                 :            :     {
    1021                 :      15470 :       GEN t = gsubgs(gmulgs(n,j),i-j);
    1022                 :      15470 :       s = gadd(s, gmul(gmul(t, gel(X,j)), gel(Y,i-j)));
    1023                 :            :     }
    1024                 :       1722 :     gel(Y,i) = gerepileupto(av, gdivgs(s,i));
    1025                 :            :   }
    1026                 :        161 :   return y;
    1027                 :            : }
    1028                 :            : 
    1029                 :            : /* we suppose n != 0, valp(x) = 0 and leading-term(x) != 0. Not stack clean */
    1030                 :            : static GEN
    1031                 :        161 : ser_pow(GEN x, GEN n, long prec)
    1032                 :            : {
    1033                 :            :   GEN y, c, lead;
    1034         [ -  + ]:        161 :   if (varncmp(gvar(n), varn(x)) <= 0) return gexp(gmul(n, glog(x,prec)), prec);
    1035                 :        161 :   lead = gel(x,2);
    1036         [ +  + ]:        161 :   if (gequal1(lead)) return ser_pow_1(x, n);
    1037                 :         14 :   x = ser_normalize(x);
    1038 [ +  - ][ +  - ]:         14 :   if (typ(n) == t_FRAC && !isinexact(lead) && ispower(lead, gel(n,2), &c))
                 [ +  + ]
    1039                 :          7 :     c = powgi(c, gel(n,1));
    1040                 :            :   else
    1041                 :          7 :     c = gpow(lead,n, prec);
    1042                 :         14 :   y = gmul(c, ser_pow_1(x, n));
    1043                 :            :   /* gpow(t_POLMOD,n) can be a t_COL [conjvec] */
    1044         [ -  + ]:         14 :   if (typ(y) != t_SER) pari_err_TYPE("gpow", y);
    1045                 :        161 :   return y;
    1046                 :            : }
    1047                 :            : 
    1048                 :            : static long
    1049                 :        175 : val_from_i(GEN E)
    1050                 :            : {
    1051         [ +  + ]:        175 :   if (is_bigint(E)) pari_err_OVERFLOW("sqrtn [valuation]");
    1052                 :        168 :   return itos(E);
    1053                 :            : }
    1054                 :            : 
    1055                 :            : /* return x^q, assume typ(x) = t_SER, typ(q) = t_INT/t_FRAC and q != 0 */
    1056                 :            : static GEN
    1057                 :        182 : ser_powfrac(GEN x, GEN q, long prec)
    1058                 :            : {
    1059                 :        182 :   GEN y, E = gmulsg(valp(x), q);
    1060                 :            :   long e;
    1061                 :            : 
    1062         [ +  + ]:        182 :   if (!signe(x))
    1063                 :            :   {
    1064         [ -  + ]:         21 :     if (gsigne(q) < 0) pari_err_INV("gpow", x);
    1065                 :         21 :     return zeroser(varn(x), val_from_i(gfloor(E)));
    1066                 :            :   }
    1067         [ +  + ]:        161 :   if (typ(E) != t_INT)
    1068                 :          7 :     pari_err_DOMAIN("sqrtn", "valuation", "!=", mkintmod(gen_0, gel(q,2)), x);
    1069                 :        154 :   e = val_from_i(E);
    1070                 :        154 :   y = leafcopy(x); setvalp(y, 0);
    1071                 :        154 :   y = ser_pow(y, q, prec);
    1072                 :        168 :   setvalp(y, e); return y;
    1073                 :            : }
    1074                 :            : 
    1075                 :            : GEN
    1076                 :    8220299 : gpow(GEN x, GEN n, long prec)
    1077                 :            : {
    1078                 :    8220299 :   long i, lx, tx, tn = typ(n);
    1079                 :            :   pari_sp av;
    1080                 :            :   GEN y;
    1081                 :            : 
    1082         [ +  + ]:    8220299 :   if (tn == t_INT) return powgi(x,n);
    1083                 :     152417 :   tx = typ(x);
    1084         [ +  + ]:     152417 :   if (is_matvec_t(tx))
    1085                 :            :   {
    1086                 :         28 :     y = cgetg_copy(x, &lx);
    1087         [ +  + ]:         70 :     for (i=1; i<lx; i++) gel(y,i) = gpow(gel(x,i),n,prec);
    1088                 :         28 :     return y;
    1089                 :            :   }
    1090                 :     152389 :   av = avma;
    1091      [ +  +  + ]:     152389 :   switch (tx)
    1092                 :            :   {
    1093                 :         14 :     case t_POL: case t_RFRAC: x = toser_i(x); /* fall through */
    1094                 :            :     case t_SER:
    1095         [ +  + ]:         77 :       if (tn == t_FRAC) return gerepileupto(av, ser_powfrac(x, n, prec));
    1096         [ +  + ]:         35 :       if (valp(x))
    1097                 :         21 :         pari_err_DOMAIN("gpow [irrational exponent]",
    1098                 :            :                         "valuation", "!=", gen_0, x);
    1099         [ +  + ]:         14 :       if (lg(x) == 2) return gerepilecopy(av, x); /* O(1) */
    1100                 :          7 :       return gerepileupto(av, ser_pow(x, n, prec));
    1101                 :            :   }
    1102         [ +  + ]:     152312 :   if (gequal0(x))
    1103                 :            :   {
    1104         [ +  + ]:         35 :     switch(tn)
    1105                 :            :     {
    1106                 :            :       case t_REAL: case t_FRAC: case t_COMPLEX: case t_QUAD:
    1107                 :         28 :         break;
    1108                 :          7 :       default: pari_err_TYPE("gpow(0,n)", n);
    1109                 :            :     }
    1110                 :         28 :     n = real_i(n);
    1111         [ +  + ]:         28 :     if (gsigne(n) <= 0) pari_err_DOMAIN("gpow(0,n)", "n", "<=", gen_0, n);
    1112         [ +  + ]:         21 :     if (!precision(x)) return gcopy(x);
    1113                 :            : 
    1114                 :         14 :     x = ground(gmulsg(gexpo(x),n));
    1115 [ +  + ][ -  + ]:         14 :     if (is_bigint(x) || uel(x,2) >= HIGHEXPOBIT)
    1116                 :          7 :       pari_err_OVERFLOW("gpow");
    1117                 :          7 :     avma = av; return real_0_bit(itos(x));
    1118                 :            :   }
    1119         [ +  + ]:     152277 :   if (tn == t_FRAC)
    1120                 :            :   {
    1121                 :      26495 :     GEN z, d = gel(n,2), a = gel(n,1);
    1122   [ +  +  +  + ]:      26495 :     switch (tx)
    1123                 :            :     {
    1124                 :            :     case t_INTMOD:
    1125                 :            :       {
    1126                 :         21 :         GEN p = gel(x,1);
    1127         [ +  + ]:         21 :         if (!BPSW_psp(p)) pari_err_PRIME("gpow",p);
    1128                 :         14 :         y = cgetg(3,t_INTMOD); gel(y,1) = icopy(p);
    1129                 :         14 :         av = avma;
    1130                 :         14 :         z = Fp_sqrtn(gel(x,2), d, p, NULL);
    1131         [ +  + ]:         14 :         if (!z) pari_err_SQRTN("gpow",x);
    1132                 :          7 :         gel(y,2) = gerepileuptoint(av, Fp_pow(z, a, p));
    1133                 :          7 :         return y;
    1134                 :            :       }
    1135                 :            : 
    1136                 :            :     case t_PADIC:
    1137                 :         14 :       z = Qp_sqrtn(x, d, NULL);
    1138         [ +  + ]:         14 :       if (!z) pari_err_SQRTN("gpow",x);
    1139                 :          7 :       return gerepileupto(av, powgi(z, a));
    1140                 :            : 
    1141                 :            :     case t_FFELT:
    1142                 :         21 :       return gerepileupto(av,FF_pow(FF_sqrtn(x,d,NULL),a));
    1143                 :            :     }
    1144                 :            :   }
    1145         [ +  + ]:     152221 :   i = (long) precision(n); if (i) prec=i;
    1146                 :     152221 :   y = gmul(n, glog(x,prec));
    1147                 :    8220187 :   return gerepileupto(av, gexp(y,prec));
    1148                 :            : }
    1149                 :            : 
    1150                 :            : GEN
    1151                 :        280 : gpowers(GEN x, long n)
    1152                 :            : {
    1153         [ -  + ]:        280 :   if (n < 0)
    1154                 :          0 :     pari_err_DOMAIN("powers", "n", "<", gen_0, stoi(n));
    1155                 :        280 :   return gen_powers(x, n, 1, (void*)x, &_sqr, &_mul, &_one);
    1156                 :            : }
    1157                 :            : 
    1158                 :            : /********************************************************************/
    1159                 :            : /**                                                                **/
    1160                 :            : /**                        RACINE CARREE                           **/
    1161                 :            : /**                                                                **/
    1162                 :            : /********************************************************************/
    1163                 :            : /* assume x unit, precp(x) = pp > 3 */
    1164                 :            : static GEN
    1165                 :        581 : sqrt_2adic(GEN x, long pp)
    1166                 :            : {
    1167 [ +  + ][ +  + ]:        581 :   GEN z = mod16(x)==(signe(x)>=0?1:15)?gen_1: utoipos(3);
    1168                 :            :   long zp;
    1169                 :            :   pari_sp av;
    1170                 :            : 
    1171         [ +  + ]:        581 :   if (pp == 4) return z;
    1172                 :        511 :   zp = 3; /* number of correct bits in z (compared to sqrt(x)) */
    1173                 :            : 
    1174                 :        511 :   av = avma;
    1175                 :            :   for(;;)
    1176                 :            :   {
    1177                 :            :     GEN mod;
    1178                 :       1561 :     zp = (zp<<1) - 1;
    1179         [ +  + ]:       1561 :     if (zp > pp) zp = pp;
    1180                 :       1561 :     mod = int2n(zp);
    1181                 :       1561 :     z = addii(z, remi2n(mulii(x, Fp_inv(z,mod)), zp));
    1182                 :       1561 :     z = shifti(z, -1); /* (z + x/z) / 2 */
    1183         [ +  + ]:       1561 :     if (pp == zp) return z;
    1184                 :            : 
    1185         [ +  - ]:       1050 :     if (zp < pp) zp--;
    1186                 :            : 
    1187         [ -  + ]:       1050 :     if (gc_needed(av,2))
    1188                 :            :     {
    1189         [ #  # ]:          0 :       if (DEBUGMEM > 1) pari_warn(warnmem,"Qp_sqrt");
    1190                 :          0 :       z = gerepileuptoint(av,z);
    1191                 :            :     }
    1192                 :       1631 :   }
    1193                 :            : }
    1194                 :            : 
    1195                 :            : /* x unit defined modulo p^e, e > 0 */
    1196                 :            : static GEN
    1197                 :       2296 : Up_sqrt(GEN x, GEN p, long e)
    1198                 :            : {
    1199                 :       2296 :   pari_sp av = avma;
    1200         [ +  + ]:       2296 :   if (equaliu(p,2))
    1201                 :            :   {
    1202         [ +  + ]:        721 :     long r = signe(x)>=0?mod8(x):8-mod8(x);
    1203         [ +  + ]:        721 :     if (e <= 3)
    1204                 :            :     {
    1205   [ +  +  +  - ]:        140 :       switch(e) {
    1206                 :          7 :       case 1: break;
    1207         [ +  + ]:        119 :       case 2: if ((r&3) == 1) break;
    1208                 :          7 :               return NULL;
    1209         [ +  + ]:         14 :       case 3: if (r == 1) break;
    1210                 :          7 :               return NULL;
    1211                 :            :       }
    1212                 :        126 :       return gen_1;
    1213                 :            :     }
    1214                 :            :     else
    1215                 :            :     {
    1216         [ -  + ]:        581 :       if (r != 1) return NULL;
    1217                 :        581 :       return gerepileuptoint(av, sqrt_2adic(x, e));
    1218                 :            :     }
    1219                 :            :   }
    1220                 :            :   else
    1221                 :       2296 :     return Zp_sqrt(x, p, e);
    1222                 :            : }
    1223                 :            : 
    1224                 :            : GEN
    1225                 :       1358 : Qp_sqrt(GEN x)
    1226                 :            : {
    1227                 :       1358 :   long pp, e = valp(x);
    1228                 :       1358 :   GEN z,y,mod, p = gel(x,2);
    1229                 :            : 
    1230         [ -  + ]:       1358 :   if (gequal0(x)) return zeropadic(p, (e+1) >> 1);
    1231         [ +  + ]:       1358 :   if (e & 1) return NULL;
    1232                 :            : 
    1233                 :       1344 :   y = cgetg(5,t_PADIC);
    1234                 :       1344 :   pp = precp(x);
    1235                 :       1344 :   mod = gel(x,3);
    1236                 :       1344 :   z   = gel(x,4); /* lift to t_INT */
    1237                 :       1344 :   e >>= 1;
    1238                 :       1344 :   z = Up_sqrt(z, p, pp);
    1239         [ +  + ]:       1344 :   if (!z) return NULL;
    1240         [ +  + ]:       1302 :   if (equaliu(p,2))
    1241                 :            :   {
    1242         [ +  + ]:        399 :     pp  = (pp <= 3) ? 1 : pp-1;
    1243                 :        399 :     mod = int2n(pp);
    1244                 :            :   }
    1245                 :        903 :   else mod = icopy(mod);
    1246                 :       1302 :   y[1] = evalprecp(pp) | evalvalp(e);
    1247                 :       1302 :   gel(y,2) = icopy(p);
    1248                 :       1302 :   gel(y,3) = mod;
    1249                 :       1358 :   gel(y,4) = z; return y;
    1250                 :            : }
    1251                 :            : 
    1252                 :            : GEN
    1253                 :        350 : Zn_sqrt(GEN d, GEN fn)
    1254                 :            : {
    1255                 :        350 :   pari_sp ltop = avma, btop;
    1256                 :        350 :   GEN b = gen_0, m = gen_1;
    1257                 :            :   long j, np;
    1258         [ -  + ]:        350 :   if (typ(d) != t_INT) pari_err_TYPE("Zn_sqrt",d);
    1259         [ -  + ]:        350 :   if (typ(fn) == t_INT)
    1260                 :          0 :     fn = absi_factor(fn);
    1261         [ -  + ]:        350 :   else if (!is_Z_factorpos(fn))
    1262                 :          0 :     pari_err_TYPE("Zn_sqrt",fn);
    1263                 :        350 :   np = nbrows(fn);
    1264                 :        350 :   btop = avma;
    1265         [ +  + ]:       1400 :   for (j = 1; j <= np; ++j)
    1266                 :            :   {
    1267                 :            :     GEN  bp, mp, pr, r;
    1268                 :       1050 :     GEN  p = gcoeff(fn, j, 1);
    1269                 :       1050 :     long e = itos(gcoeff(fn, j, 2));
    1270                 :       1050 :     long v = Z_pvalrem(d,p,&r);
    1271         [ +  + ]:       1050 :     if (v >= e) bp =gen_0;
    1272                 :            :     else
    1273                 :            :     {
    1274         [ -  + ]:        952 :       if (odd(v)) return NULL;
    1275                 :        952 :       bp = Up_sqrt(r, p, e-v);
    1276         [ -  + ]:        952 :       if (!bp)    return NULL;
    1277         [ -  + ]:        952 :       if (v) bp = mulii(bp, powiu(p, v>>1L));
    1278                 :            :     }
    1279                 :       1050 :     mp = powiu(p, e);
    1280                 :       1050 :     pr = mulii(m, mp);
    1281                 :       1050 :     b = Z_chinese_coprime(b, bp, m, mp, pr);
    1282                 :       1050 :     m = pr;
    1283         [ -  + ]:       1050 :     if (gc_needed(btop, 1))
    1284                 :          0 :       gerepileall(btop, 2, &b, &m);
    1285                 :            :   }
    1286                 :        350 :   return gerepileupto(ltop, b);
    1287                 :            : }
    1288                 :            : 
    1289                 :            : static GEN
    1290                 :        973 : sqrt_ser(GEN b, long prec)
    1291                 :            : {
    1292                 :        973 :   long e = valp(b), vx = varn(b), lx, lold, j;
    1293                 :            :   ulong mask;
    1294                 :            :   GEN a, x, lta, ltx;
    1295                 :            : 
    1296         [ -  + ]:        973 :   if (!signe(b)) return zeroser(vx, e>>1);
    1297                 :        973 :   a = leafcopy(b);
    1298                 :        973 :   x = cgetg_copy(b, &lx);
    1299         [ +  + ]:        973 :   if (e & 1)
    1300                 :          7 :     pari_err_DOMAIN("sqrtn", "valuation", "!=", mkintmod(gen_0, gen_2), x);
    1301                 :        966 :   a[1] = x[1] = evalsigne(1) | evalvarn(0) | _evalvalp(0);
    1302                 :        966 :   lta = gel(a,2);
    1303         [ +  + ]:        966 :   if (gequal1(lta)) ltx = lta;
    1304         [ +  + ]:         28 :   else if (!issquareall(lta,&ltx)) ltx = gsqrt(lta,prec);
    1305                 :        966 :   gel(x,2) = ltx;
    1306         [ +  + ]:       7154 :   for (j = 3; j < lx; j++) gel(x,j) = gen_0;
    1307                 :        966 :   setlg(x,3);
    1308                 :        966 :   mask = quadratic_prec_mask(lx - 2);
    1309                 :        966 :   lold = 1;
    1310         [ +  + ]:       3402 :   while (mask > 1)
    1311                 :            :   {
    1312                 :       2436 :     GEN y, x2 = gmul2n(x,1);
    1313                 :       2436 :     long l = lold << 1;
    1314                 :            : 
    1315         [ +  + ]:       2436 :     if (mask & 1) l--;
    1316                 :       2436 :     mask >>= 1;
    1317                 :       2436 :     setlg(a, l + 2);
    1318                 :       2436 :     setlg(x, l + 2);
    1319                 :       2436 :     y = sqr_ser_part(x, lold, l-1) - lold;
    1320         [ +  + ]:       8624 :     for (j = lold+2; j < l+2; j++) gel(y,j) = gsub(gel(y,j), gel(a,j));
    1321                 :       2436 :     y += lold; setvalp(y, lold);
    1322                 :       2436 :     y = normalize(y);
    1323                 :       2436 :     y = gsub(x, gdiv(y, x2)); /* = gmul2n(gsub(x, gdiv(a,x)), -1); */
    1324         [ +  + ]:       8624 :     for (j = lold+2; j < l+2; j++) gel(x,j) = gel(y,j);
    1325                 :       2436 :     lold = l;
    1326                 :            :   }
    1327                 :        966 :   x[1] = evalsigne(1) | evalvarn(vx) | _evalvalp(e >> 1);
    1328                 :        966 :   return x;
    1329                 :            : }
    1330                 :            : 
    1331                 :            : GEN
    1332                 :    3358333 : gsqrt(GEN x, long prec)
    1333                 :            : {
    1334                 :            :   pari_sp av;
    1335                 :            :   GEN y;
    1336                 :            : 
    1337   [ +  +  +  +  :    3358333 :   switch(typ(x))
                   +  + ]
    1338                 :            :   {
    1339                 :    2650808 :     case t_REAL: return sqrtr(x);
    1340                 :            : 
    1341                 :            :     case t_INTMOD:
    1342                 :            :     {
    1343                 :         35 :       GEN p = gel(x,1), a;
    1344                 :         35 :       y = cgetg(3,t_INTMOD); gel(y,1) = icopy(p);
    1345                 :         35 :       a = Fp_sqrt(gel(x,2),p);
    1346         [ +  + ]:         21 :       if (!a)
    1347                 :            :       {
    1348         [ -  + ]:          7 :         if (!BPSW_psp(p)) pari_err_PRIME("sqrt [modulus]",p);
    1349                 :          7 :         pari_err_SQRTN("gsqrt",x);
    1350                 :            :       }
    1351                 :         14 :       gel(y,2) = a; return y;
    1352                 :            :     }
    1353                 :            : 
    1354                 :            :     case t_COMPLEX:
    1355                 :            :     { /* (u+iv)^2 = a+ib <=> u^2+v^2 = sqrt(a^2+b^2), u^2-v^2=a, 2uv=b */
    1356                 :     440226 :       GEN a = gel(x,1), b = gel(x,2), r, u, v;
    1357         [ -  + ]:     440226 :       if (isrationalzero(b)) return gsqrt(a, prec);
    1358                 :     440226 :       y = cgetg(3,t_COMPLEX); av = avma;
    1359                 :            : 
    1360                 :     440226 :       r = cxnorm(x);
    1361         [ -  + ]:     440226 :       if (typ(r) == t_INTMOD) pari_err_IMPL("sqrt(complex of t_INTMODs)");
    1362                 :     440226 :       r = gsqrt(r, prec); /* t_REAL, |a+Ib| */
    1363         [ +  + ]:     440226 :       if (!signe(r))
    1364                 :         56 :         u = v = gerepileuptoleaf(av, sqrtr(r));
    1365         [ +  + ]:     440170 :       else if (gsigne(a) < 0)
    1366                 :            :       {
    1367                 :            :         /* v > 0 since r > 0, a < 0, rounding errors can't make the sum of two
    1368                 :            :          * positive numbers = 0 */
    1369                 :      39092 :         v = sqrtr( gmul2n(gsub(r,a), -1) );
    1370         [ +  + ]:      39092 :         if (gsigne(b) < 0) togglesign(v);
    1371                 :      39092 :         v = gerepileuptoleaf(av, v); av = avma;
    1372                 :            :         /* v = 0 is impossible */
    1373                 :      39092 :         u = gerepileuptoleaf(av, gdiv(b, shiftr(v,1)));
    1374                 :            :       } else {
    1375                 :     401078 :         u = sqrtr( gmul2n(gadd(r,a), -1) );
    1376                 :     401078 :         u = gerepileuptoleaf(av, u); av = avma;
    1377         [ +  + ]:     401078 :         if (!signe(u)) /* possible if a = 0.0, e.g. sqrt(0.e-10+1e-10*I) */
    1378                 :          7 :           v = u;
    1379                 :            :         else
    1380                 :     401071 :           v = gerepileuptoleaf(av, gdiv(b, shiftr(u,1)));
    1381                 :            :       }
    1382                 :     440226 :       gel(y,1) = u;
    1383                 :     440226 :       gel(y,2) = v; return y;
    1384                 :            :     }
    1385                 :            : 
    1386                 :            :     case t_PADIC:
    1387                 :         49 :       y = Qp_sqrt(x);
    1388         [ +  + ]:         49 :       if (!y) pari_err_SQRTN("Qp_sqrt",x);
    1389                 :         28 :       return y;
    1390                 :            : 
    1391                 :         70 :     case t_FFELT: return FF_sqrt(x);
    1392                 :            : 
    1393                 :            :     default:
    1394         [ +  + ]:     267145 :       av = avma; if (!(y = toser_i(x))) break;
    1395                 :        973 :       return gerepilecopy(av, sqrt_ser(y, prec));
    1396                 :            :   }
    1397                 :    3358284 :   return trans_eval("sqrt",gsqrt,x,prec);
    1398                 :            : }
    1399                 :            : /********************************************************************/
    1400                 :            : /**                                                                **/
    1401                 :            : /**                          N-th ROOT                             **/
    1402                 :            : /**                                                                **/
    1403                 :            : /********************************************************************/
    1404                 :            : /* exp(2Ipi/n), assume n positive t_INT */
    1405                 :            : GEN
    1406                 :          7 : rootsof1complex(GEN n, long prec)
    1407                 :            : {
    1408                 :          7 :   pari_sp av = avma;
    1409         [ -  + ]:          7 :   if (is_pm1(n)) return real_1(prec);
    1410         [ -  + ]:          7 :   if (equaliu(n, 2)) return stor(-1, prec);
    1411                 :          7 :   return gerepileupto(av, expIr( divri(Pi2n(1, prec), n) ));
    1412                 :            : }
    1413                 :            : 
    1414                 :            : /*Only the O() of y is used*/
    1415                 :            : GEN
    1416                 :          0 : rootsof1padic(GEN n, GEN y)
    1417                 :            : {
    1418                 :          0 :   GEN z, r = cgetp(y), p = gel(y,2);
    1419                 :          0 :   pari_sp av = avma;
    1420                 :          0 :   z = rootsof1_Fp(n, p);
    1421                 :          0 :   z = Zp_sqrtnlift(gen_1, n, z, p, precp(y));
    1422                 :          0 :   affii(z, gel(r,4)); avma = av; return r;
    1423                 :            : }
    1424                 :            : 
    1425                 :            : /* Let x = 1 mod p and y := (x-1)/(x+1) = 0 (p). Then
    1426                 :            :  * log(x) = log(1+y) - log(1-y) = 2 \sum_{k odd} y^k / k.
    1427                 :            :  * palogaux(x) returns the last sum (not multiplied by 2) */
    1428                 :            : static GEN
    1429                 :       1057 : palogaux(GEN x)
    1430                 :            : {
    1431                 :            :   long k,e,pp;
    1432                 :       1057 :   GEN y,s,y2, p = gel(x,2);
    1433                 :            : 
    1434         [ +  + ]:       1057 :   if (equalii(gen_1, gel(x,4)))
    1435                 :            :   {
    1436                 :        420 :     long v = valp(x)+precp(x);
    1437         [ +  + ]:        420 :     if (equaliu(p,2)) v--;
    1438                 :        420 :     return zeropadic(p, v);
    1439                 :            :   }
    1440                 :        637 :   y = gdiv(gaddgs(x,-1), gaddgs(x,1));
    1441                 :        637 :   e = valp(y); /* > 0 */
    1442         [ +  + ]:        637 :   if (e <= 0) {
    1443         [ +  - ]:          7 :     if (!BPSW_psp(p)) pari_err_PRIME("p-adic log",p);
    1444                 :          0 :     pari_err_BUG("log_p");
    1445                 :            :   }
    1446                 :        630 :   pp = e+precp(y);
    1447         [ +  + ]:        630 :   if (equaliu(p,2)) pp--;
    1448                 :            :   else
    1449                 :            :   {
    1450                 :            :     GEN p1;
    1451         [ +  + ]:       1309 :     for (p1=utoipos(e); cmpui(pp,p1) > 0; pp++) p1 = mulii(p1, p);
    1452                 :        588 :     pp -= 2;
    1453                 :            :   }
    1454         [ +  + ]:        630 :   k = pp/e; if (!odd(k)) k--;
    1455                 :        630 :   y2 = gsqr(y); s = gdivgs(gen_1,k);
    1456         [ +  + ]:       1001 :   while (k > 2)
    1457                 :            :   {
    1458                 :        371 :     k -= 2; s = gadd(gmul(y2,s), gdivgs(gen_1,k));
    1459                 :            :   }
    1460                 :       1050 :   return gmul(s,y);
    1461                 :            : }
    1462                 :            : 
    1463                 :            : GEN
    1464                 :       1064 : Qp_log(GEN x)
    1465                 :            : {
    1466                 :       1064 :   pari_sp av = avma;
    1467                 :       1064 :   GEN y, p = gel(x,2), a = gel(x,4);
    1468                 :            : 
    1469         [ +  + ]:       1064 :   if (!signe(a)) pari_err_DOMAIN("Qp_log", "argument", "=", gen_0, x);
    1470                 :       1057 :   y = leafcopy(x); setvalp(y,0);
    1471         [ +  + ]:       1057 :   if (equaliu(p,2))
    1472                 :        280 :     y = palogaux(gsqr(y));
    1473         [ +  + ]:        777 :   else if (gequal1(modii(a, p)))
    1474                 :        287 :     y = gmul2n(palogaux(y), 1);
    1475                 :            :   else
    1476                 :            :   { /* compute log(x^(p-1)) / (p-1) */
    1477                 :        490 :     GEN mod = gel(y,3), p1 = subis(p,1);
    1478                 :        490 :     gel(y,4) = Fp_pow(a, p1, mod);
    1479                 :        490 :     p1 = diviiexact(subsi(1,mod), p1); /* 1/(p-1) */
    1480                 :        490 :     y = gmul(palogaux(y), shifti(p1,1));
    1481                 :            :   }
    1482                 :       1050 :   return gerepileupto(av,y);
    1483                 :            : }
    1484                 :            : 
    1485                 :            : static GEN Qp_exp_safe(GEN x);
    1486                 :            : 
    1487                 :            : /*compute the p^e th root of x p-adic, assume x != 0 */
    1488                 :            : static GEN
    1489                 :        854 : Qp_sqrtn_ram(GEN x, long e)
    1490                 :            : {
    1491                 :        854 :   pari_sp ltop=avma;
    1492                 :        854 :   GEN a, p = gel(x,2), n = powiu(p,e);
    1493                 :        854 :   long v = valp(x), va;
    1494         [ +  + ]:        854 :   if (v)
    1495                 :            :   {
    1496                 :            :     long z;
    1497                 :        161 :     v = sdivsi_rem(v, n, &z);
    1498         [ +  + ]:        161 :     if (z) return NULL;
    1499                 :         91 :     x = leafcopy(x);
    1500                 :         91 :     setvalp(x,0);
    1501                 :            :   }
    1502                 :            :   /*If p = 2, -1 is a root of 1 in U1: need extra check*/
    1503 [ +  + ][ +  + ]:        784 :   if (equaliu(p, 2) && mod8(gel(x,4)) != 1) return NULL;
    1504                 :        749 :   a = Qp_log(x);
    1505                 :        749 :   va = valp(a) - e;
    1506         [ +  + ]:        749 :   if (va <= 0)
    1507                 :            :   {
    1508         [ +  + ]:        287 :     if (signe(gel(a,4))) return NULL;
    1509                 :            :     /* all accuracy lost */
    1510                 :        119 :     a = cvtop(remii(gel(x,4),p), p, 1);
    1511                 :            :   }
    1512                 :            :   else
    1513                 :            :   {
    1514                 :        462 :     setvalp(a, va); /* divide by p^e */
    1515                 :        462 :     a = Qp_exp_safe(a);
    1516         [ -  + ]:        462 :     if (!a) return NULL;
    1517                 :            :     /* n=p^e and a^n=z*x where z is a (p-1)th-root of 1.
    1518                 :            :      * Since z^n=z, we have (a/z)^n = x. */
    1519                 :        462 :     a = gdiv(x, powgi(a,addis(n,-1))); /* = a/z = x/a^(n-1)*/
    1520         [ +  + ]:        462 :     if (v) setvalp(a,v);
    1521                 :            :   }
    1522                 :        854 :   return gerepileupto(ltop,a);
    1523                 :            : }
    1524                 :            : 
    1525                 :            : /*compute the nth root of x p-adic p prime with n*/
    1526                 :            : static GEN
    1527                 :        616 : Qp_sqrtn_unram(GEN x, GEN n, GEN *zetan)
    1528                 :            : {
    1529                 :            :   pari_sp av;
    1530                 :        616 :   GEN Z, a, r, p = gel(x,2);
    1531                 :        616 :   long v = valp(x);
    1532         [ +  + ]:        616 :   if (v)
    1533                 :            :   {
    1534                 :            :     long z;
    1535                 :         84 :     v = sdivsi_rem(v,n,&z);
    1536         [ +  + ]:         84 :     if (z) return NULL;
    1537                 :            :   }
    1538                 :        609 :   r = cgetp(x); setvalp(r,v);
    1539                 :        609 :   Z = NULL; /* -Wall */
    1540         [ +  + ]:        609 :   if (zetan) Z = cgetp(x);
    1541                 :        609 :   av = avma; a = Fp_sqrtn(gel(x,4), n, p, zetan);
    1542         [ +  + ]:        609 :   if (!a) return NULL;
    1543                 :        595 :   affii(Zp_sqrtnlift(gel(x,4), n, a, p, precp(x)), gel(r,4));
    1544         [ +  + ]:        595 :   if (zetan)
    1545                 :            :   {
    1546                 :         14 :     affii(Zp_sqrtnlift(gen_1, n, *zetan, p, precp(x)), gel(Z,4));
    1547                 :         14 :     *zetan = Z;
    1548                 :            :   }
    1549                 :        616 :   avma = av; return r;
    1550                 :            : }
    1551                 :            : 
    1552                 :            : GEN
    1553                 :       1183 : Qp_sqrtn(GEN x, GEN n, GEN *zetan)
    1554                 :            : {
    1555                 :            :   pari_sp av, tetpil;
    1556                 :            :   GEN q, p;
    1557                 :            :   long e;
    1558         [ +  + ]:       1183 :   if (equaliu(n, 2))
    1559                 :            :   {
    1560         [ -  + ]:         70 :     if (zetan) *zetan = gen_m1;
    1561         [ +  + ]:         70 :     if (signe(n) < 0) x = ginv(x);
    1562                 :         63 :     return Qp_sqrt(x);
    1563                 :            :   }
    1564                 :       1113 :   av = avma; p = gel(x,2);
    1565         [ +  + ]:       1113 :   if (!signe(gel(x,4)))
    1566                 :            :   {
    1567         [ -  + ]:        203 :     if (signe(n) < 0) pari_err_INV("Qp_sqrtn", x);
    1568                 :        203 :     q = divii(addis(n, valp(x)-1), n);
    1569         [ -  + ]:        203 :     if (zetan) *zetan = gen_1;
    1570                 :        203 :     avma = av; return zeropadic(p, itos(q));
    1571                 :            :   }
    1572                 :            :   /* treat the ramified part using logarithms */
    1573                 :        910 :   e = Z_pvalrem(n, p, &q);
    1574 [ +  + ][ +  + ]:        910 :   if (e) { x = Qp_sqrtn_ram(x,e); if (!x) return NULL; }
    1575         [ +  + ]:        637 :   if (is_pm1(q))
    1576                 :            :   { /* finished */
    1577         [ +  + ]:         21 :     if (signe(q) < 0) x = ginv(x);
    1578                 :         21 :     x = gerepileupto(av, x);
    1579         [ +  + ]:         21 :     if (zetan)
    1580         [ +  + ]:         28 :       *zetan = (e && equaliu(p, 2))? gen_m1 /*-1 in Q_2*/
    1581         [ +  - ]:         28 :                                    : gen_1;
    1582                 :         21 :     return x;
    1583                 :            :   }
    1584                 :        616 :   tetpil = avma;
    1585                 :            :   /* use hensel lift for unramified case */
    1586                 :        616 :   x = Qp_sqrtn_unram(x, q, zetan);
    1587         [ +  + ]:        616 :   if (!x) return NULL;
    1588         [ +  + ]:        595 :   if (zetan)
    1589                 :            :   {
    1590                 :            :     GEN *gptr[2];
    1591 [ +  + ][ +  - ]:         14 :     if (e && equaliu(p, 2))/*-1 in Q_2*/
    1592                 :            :     {
    1593                 :          7 :       tetpil = avma; x = gcopy(x); *zetan = gneg(*zetan);
    1594                 :            :     }
    1595                 :         14 :     gptr[0] = &x; gptr[1] = zetan;
    1596                 :         14 :     gerepilemanysp(av,tetpil,gptr,2);
    1597                 :         14 :     return x;
    1598                 :            :   }
    1599                 :       1176 :   return gerepile(av,tetpil,x);
    1600                 :            : }
    1601                 :            : 
    1602                 :            : GEN
    1603                 :       1680 : sqrtnint(GEN a, long n)
    1604                 :            : {
    1605                 :       1680 :   pari_sp ltop = avma;
    1606                 :            :   GEN x, b, q;
    1607                 :            :   long s, k, e;
    1608                 :       1680 :   const ulong nm1 = n - 1;
    1609         [ -  + ]:       1680 :   if (typ(a) != t_INT) pari_err_TYPE("sqrtnint",a);
    1610         [ +  + ]:       1680 :   if (n <= 0) pari_err_DOMAIN("sqrtnint", "n", "<=", gen_0, stoi(n));
    1611         [ -  + ]:       1673 :   if (n == 1) return icopy(a);
    1612                 :       1673 :   s = signe(a);
    1613         [ +  + ]:       1673 :   if (s < 0) pari_err_DOMAIN("sqrtnint", "x", "<", gen_0, a);
    1614         [ -  + ]:       1666 :   if (!s) return gen_0;
    1615         [ +  + ]:       1666 :   if (lgefint(a) == 3) return utoi(usqrtn(itou(a), n));
    1616                 :       1237 :   e = expi(a); k = e/(2*n);
    1617         [ +  + ]:       1237 :   if (k == 0)
    1618                 :            :   {
    1619                 :            :     long flag;
    1620         [ -  + ]:        291 :     if (n > e) {avma = ltop; return gen_1;}
    1621                 :        291 :     flag = cmpii(a, powuu(3, n)); avma = ltop;
    1622         [ -  + ]:        291 :     return (flag < 0) ? gen_2: stoi(3);
    1623                 :            :   }
    1624         [ +  + ]:        946 :   if (e < n*(BITS_IN_LONG - 1))
    1625                 :            :   {
    1626                 :            :     ulong s, xs, qs;
    1627                 :        799 :     s = 1 + e/n; xs = 1UL << s;
    1628                 :        799 :     qs = itou(shifti(a, -nm1*s));
    1629         [ +  + ]:       4845 :     while (qs < xs) {
    1630                 :       4053 :       xs -= (xs - qs + nm1)/n;
    1631                 :       4053 :       q = divii(a, powuu(xs, nm1));
    1632         [ +  + ]:       4053 :       if (lgefint(q) > 3) break;
    1633                 :       4046 :       qs = itou(q);
    1634                 :            :     }
    1635                 :        799 :     return utoi(xs);
    1636                 :            :   }
    1637                 :        147 :   b = addsi(1, shifti(a, -n*k));
    1638                 :        147 :   x = shifti(addsi(1, sqrtnint(b, n)), k);
    1639                 :        147 :   q = divii(a, powiu(x, nm1));
    1640         [ +  + ]:        372 :   while (cmpii(q, x) < 0) /* a priori one iteration, no GC necessary */
    1641                 :            :   {
    1642                 :        225 :     x = subii(x, divis(addsi(nm1, subii(x, q)), n));
    1643                 :        225 :     q = divii(a, powiu(x, nm1));
    1644                 :            :   }
    1645                 :       1666 :   return gerepileuptoleaf(ltop, x);
    1646                 :            : }
    1647                 :            : 
    1648                 :            : ulong
    1649                 :        429 : usqrtn(ulong a, ulong n)
    1650                 :            : {
    1651                 :            :   ulong x, s, q;
    1652                 :        429 :   const ulong nm1 = n - 1;
    1653         [ -  + ]:        429 :   if (!n) pari_err_DOMAIN("sqrtnint", "n", "=", gen_0, utoi(n));
    1654 [ +  - ][ -  + ]:        429 :   if (n == 1 || a == 0) return a;
    1655                 :        429 :   s = 1 + expu(a)/n; x = 1UL << s;
    1656         [ +  + ]:        429 :   q = (nm1*s >= BITS_IN_LONG)? 0: a >> (nm1*s);
    1657         [ +  + ]:       1265 :   while (q < x) {
    1658                 :            :     ulong X;
    1659                 :        836 :     x -= (x - q + nm1)/n;
    1660                 :        836 :     X = upowuu(x, nm1);
    1661         [ +  + ]:        836 :     q = X? a/X: 0;
    1662                 :            :   }
    1663                 :        429 :   return x;
    1664                 :            : }
    1665                 :            : 
    1666                 :            : GEN
    1667                 :       2261 : gsqrtn(GEN x, GEN n, GEN *zetan, long prec)
    1668                 :            : {
    1669                 :            :   long i, lx, tx;
    1670                 :            :   pari_sp av;
    1671                 :            :   GEN y, z;
    1672         [ -  + ]:       2261 :   if (typ(n)!=t_INT) pari_err_TYPE("sqrtn",n);
    1673         [ -  + ]:       2261 :   if (!signe(n)) pari_err_DOMAIN("sqrtn", "n", "=", gen_0, n);
    1674         [ +  + ]:       2261 :   if (is_pm1(n))
    1675                 :            :   {
    1676         [ +  + ]:         14 :     if (zetan) *zetan = gen_1;
    1677         [ +  + ]:         14 :     return (signe(n) > 0)? gcopy(x): ginv(x);
    1678                 :            :   }
    1679         [ +  + ]:       2247 :   if (zetan) *zetan = gen_0;
    1680                 :       2247 :   tx = typ(x);
    1681         [ +  + ]:       2247 :   if (is_matvec_t(tx))
    1682                 :            :   {
    1683                 :          7 :     y = cgetg_copy(x, &lx);
    1684         [ +  + ]:         21 :     for (i=1; i<lx; i++) gel(y,i) = gsqrtn(gel(x,i),n,NULL,prec);
    1685                 :          7 :     return y;
    1686                 :            :   }
    1687                 :       2240 :   av = avma;
    1688   [ +  +  +  +  :       2240 :   switch(tx)
                   +  + ]
    1689                 :            :   {
    1690                 :            :   case t_INTMOD:
    1691                 :            :     {
    1692                 :         56 :       GEN p = gel(x,1), s;
    1693                 :         56 :       z = gen_0;
    1694                 :         56 :       y = cgetg(3,t_INTMOD);  gel(y,1) = icopy(p);
    1695         [ +  + ]:         56 :       if (zetan) { z = cgetg(3,t_INTMOD); gel(z,1) = gel(y,1); }
    1696                 :         56 :       s = Fp_sqrtn(gel(x,2),n,p,zetan);
    1697         [ +  + ]:         35 :       if (!s) {
    1698         [ +  + ]:         21 :         if (zetan) {avma=av; return gen_0;}
    1699         [ +  + ]:         14 :         if (!BPSW_psp(p)) pari_err_PRIME("sqrtn [modulus]",p);
    1700                 :          7 :         pari_err_SQRTN("gsqrtn",x);
    1701                 :            :       }
    1702                 :         14 :       gel(y,2) = s;
    1703         [ +  + ]:         14 :       if (zetan) { gel(z,2) = *zetan; *zetan = z; }
    1704                 :         14 :       return y;
    1705                 :            :     }
    1706                 :            : 
    1707                 :            :   case t_PADIC:
    1708                 :         56 :     y = Qp_sqrtn(x,n,zetan);
    1709         [ +  + ]:         49 :     if (!y) {
    1710         [ -  + ]:          7 :       if (zetan) return gen_0;
    1711                 :          7 :       pari_err_SQRTN("gsqrtn",x);
    1712                 :            :     }
    1713                 :         42 :     return y;
    1714                 :            : 
    1715                 :         77 :   case t_FFELT: return FF_sqrtn(x,n,zetan);
    1716                 :            : 
    1717                 :            :   case t_INT: case t_FRAC: case t_REAL: case t_COMPLEX:
    1718         [ +  + ]:       1904 :     i = precision(x); if (i) prec = i;
    1719         [ +  + ]:       1904 :     if (isint1(x))
    1720                 :          7 :       y = real_1(prec);
    1721         [ +  + ]:       1897 :     else if (gequal0(x))
    1722                 :            :     {
    1723                 :            :       long b;
    1724         [ -  + ]:         21 :       if (signe(n) < 0) pari_err_INV("gsqrtn",x);
    1725         [ +  + ]:         21 :       if (isinexactreal(x))
    1726                 :         14 :         b = sdivsi(gexpo(x), n);
    1727                 :            :       else
    1728                 :          7 :         b = -prec2nbits(prec);
    1729         [ +  + ]:         21 :       if (typ(x) == t_COMPLEX)
    1730                 :            :       {
    1731                 :          7 :         y = cgetg(3,t_COMPLEX);
    1732                 :          7 :         gel(y,1) = gel(y,2) = real_0_bit(b);
    1733                 :            :       }
    1734                 :            :       else
    1735                 :         14 :         y = real_0_bit(b);
    1736                 :            :     }
    1737                 :            :     else
    1738                 :       1876 :       y = gerepileupto(av, gexp(gdiv(glog(x,prec), n), prec));
    1739         [ +  + ]:       1904 :     if (zetan) *zetan = rootsof1complex(n,prec);
    1740                 :       1904 :     return y;
    1741                 :            : 
    1742                 :            :   case t_QUAD:
    1743                 :          7 :     return gsqrtn(quadtofp(x, prec), n, zetan, prec);
    1744                 :            : 
    1745                 :            :   default:
    1746         [ -  + ]:        140 :     av = avma; if (!(y = toser_i(x))) break;
    1747                 :        140 :     return gerepileupto(av, ser_powfrac(y, ginv(n), prec));
    1748                 :            :   }
    1749                 :          0 :   pari_err_TYPE("sqrtn",x);
    1750                 :       2184 :   return NULL;/* not reached */
    1751                 :            : }
    1752                 :            : 
    1753                 :            : /********************************************************************/
    1754                 :            : /**                                                                **/
    1755                 :            : /**                             EXP(X) - 1                         **/
    1756                 :            : /**                                                                **/
    1757                 :            : /********************************************************************/
    1758                 :            : /* exp(|x|) - 1, assume x != 0.
    1759                 :            :  * For efficiency, x should be reduced mod log(2): if so, we have a < 0 */
    1760                 :            : GEN
    1761                 :    4770625 : exp1r_abs(GEN x)
    1762                 :            : {
    1763                 :    4770625 :   long l = realprec(x), a = expo(x), b = prec2nbits(l), L, i, n, m, B;
    1764                 :            :   GEN y, p2, X;
    1765                 :            :   pari_sp av;
    1766                 :            :   double d;
    1767                 :            : 
    1768         [ +  + ]:    4770625 :   if (b + a <= 0) return mpabs(x);
    1769                 :            : 
    1770                 :    4765070 :   y = cgetr(l); av = avma;
    1771                 :    4765070 :   B = b/3 + BITS_IN_LONG + (BITS_IN_LONG*BITS_IN_LONG)/ b;
    1772                 :    4765070 :   d = a/2.; m = (long)(d + sqrt(d*d + B)); /* >= 0 */
    1773         [ +  + ]:    4765070 :   if (m < (-a) * 0.1) m = 0; /* not worth it */
    1774                 :    4765070 :   L = l + nbits2extraprec(m);
    1775                 :            :  /* Multiplication is quadratic in this range (l is small, otherwise we
    1776                 :            :   * use logAGM + Newton). Set Y = 2^(-e-a) x, compute truncated series
    1777                 :            :   * sum x^k/k!: this costs roughly
    1778                 :            :   *    m b^2 + sum_{k <= n} (k e + BITS_IN_LONG)^2
    1779                 :            :   * bit operations with |x| <  2^(1+a), |Y| < 2^(1-e) , m = e+a and b bits of
    1780                 :            :   * accuracy needed, so
    1781                 :            :   *    B := (b / 3 + BITS_IN_LONG + BITS_IN_LONG^2 / b) ~ m(m-a)
    1782                 :            :   * we want b ~ 3 m (m-a) or m~b+a hence
    1783                 :            :   *     m = min( a/2 + sqrt(a^2/4 + B),  b + a )
    1784                 :            :   * NB: e ~ (b/3)^(1/2) as b -> oo
    1785                 :            :   *
    1786                 :            :   * Truncate the sum at k = n (>= 1), the remainder is
    1787                 :            :   *   sum_{k >= n+1} Y^k / k! < Y^(n+1) / (n+1)! (1-Y) < Y^(n+1) / n!
    1788                 :            :   * We want Y^(n+1) / n! <= Y 2^-b, hence -n log_2 |Y| + log_2 n! >= b
    1789                 :            :   *   log n! ~ (n + 1/2) log(n+1) - (n+1) + log(2Pi)/2,
    1790                 :            :   * error bounded by 1/6(n+1) <= 1/12. Finally, we want
    1791                 :            :   * n (-1/log(2) -log_2 |Y| + log_2(n+1)) >= b  */
    1792                 :    4765070 :   b += m;
    1793                 :    4765070 :   d = m-dbllog2(x)-1/LOG2; /* ~ -log_2 Y - 1/log(2) */
    1794                 :    4765070 :   n = (long)(b / d);
    1795         [ +  + ]:    4765070 :   if (n > 1)
    1796                 :    4753058 :     n = (long)(b / (d + log2((double)n+1))); /* log~constant in small ranges */
    1797         [ +  + ]:   10929722 :   while (n*(d+log2((double)n+1)) < b) n++; /* expect few corrections */
    1798                 :            : 
    1799                 :    4765070 :   X = rtor(x,L); shiftr_inplace(X, -m); setsigne(X, 1);
    1800         [ -  + ]:    4765070 :   if (n == 1) p2 = X;
    1801                 :            :   else
    1802                 :            :   {
    1803                 :    4765070 :     long s = 0, l1 = nbits2prec((long)(d + n + 16));
    1804                 :    4765070 :     GEN unr = real_1(L);
    1805                 :            :     pari_sp av2;
    1806                 :            : 
    1807                 :    4765070 :     p2 = cgetr(L); av2 = avma;
    1808         [ +  + ]:   78366065 :     for (i=n; i>=2; i--, avma = av2)
    1809                 :            :     { /* compute X^(n-1)/n! + ... + X/2 + 1 */
    1810                 :            :       GEN p1, p3;
    1811                 :   73600995 :       setprec(X,l1); p3 = divru(X,i);
    1812         [ +  + ]:   73600995 :       l1 += dvmdsBIL(s - expo(p3), &s); if (l1>L) l1=L;
    1813         [ +  + ]:   73600995 :       setprec(unr,l1); p1 = addrr_sign(unr,1, i == n? p3: mulrr(p3,p2),1);
    1814                 :   73600995 :       setprec(p2,l1); affrr(p1,p2); /* p2 <- 1 + (X/i)*p2 */
    1815                 :            :     }
    1816                 :    4765070 :     setprec(X,L); p2 = mulrr(X,p2);
    1817                 :            :   }
    1818                 :            : 
    1819         [ +  + ]:   55313957 :   for (i=1; i<=m; i++)
    1820                 :            :   {
    1821         [ -  + ]:   50548887 :     if (realprec(p2) > L) setprec(p2,L);
    1822                 :   50548887 :     p2 = mulrr(p2, addsr(2,p2));
    1823                 :            :   }
    1824                 :    4770625 :   affrr_fixlg(p2,y); avma = av; return y;
    1825                 :            : }
    1826                 :            : 
    1827                 :            : GEN
    1828                 :        161 : mpexpm1(GEN x)
    1829                 :            : {
    1830                 :        161 :   long sx = signe(x);
    1831                 :            :   GEN y, z;
    1832                 :            :   pari_sp av;
    1833         [ -  + ]:        161 :   if (!sx) return real_0_bit(expo(x));
    1834         [ +  + ]:        161 :   if (sx > 0) return exp1r_abs(x);
    1835                 :            :   /* compute exp(x) * (1 - exp(-x)) */
    1836                 :         21 :   av = avma; y = exp1r_abs(x);
    1837                 :         21 :   z = addsr(1, y); setsigne(z, -1);
    1838                 :        161 :   return gerepileupto(av, divrr(y, z));
    1839                 :            : }
    1840                 :            : 
    1841                 :            : GEN
    1842                 :        112 : gexpm1(GEN x, long prec)
    1843                 :            : {
    1844      [ +  +  - ]:        112 :   switch(typ(x))
    1845                 :            :   {
    1846                 :         56 :     case t_REAL: return mpexpm1(x);
    1847                 :         56 :     case t_COMPLEX: return cxexpm1(x,prec);
    1848                 :            :   }
    1849                 :        112 :   return trans_eval("expm1",gexpm1,x,prec);
    1850                 :            : }
    1851                 :            : /********************************************************************/
    1852                 :            : /**                                                                **/
    1853                 :            : /**                             EXP(X)                             **/
    1854                 :            : /**                                                                **/
    1855                 :            : /********************************************************************/
    1856                 :            : 
    1857                 :            : /* centermod(x, log(2)), set *sh to the quotient */
    1858                 :            : static GEN
    1859                 :    4685939 : modlog2(GEN x, long *sh)
    1860                 :            : {
    1861                 :    4685939 :   double d = rtodbl(x);
    1862                 :    4685939 :   long q = (long) ((fabs(d) + (LOG2/2))/LOG2);
    1863         [ -  + ]:    4685939 :   if (d > LOG2 * LONG_MAX)
    1864                 :          0 :     pari_err_OVERFLOW("expo()"); /* avoid overflow in  q */
    1865         [ +  + ]:    4685939 :   if (d < 0) q = -q;
    1866                 :    4685939 :   *sh = q;
    1867         [ +  + ]:    4685939 :   if (q) {
    1868                 :    4556861 :     long l = realprec(x) + 1;
    1869                 :    4556861 :     x = subrr(rtor(x,l), mulsr(q, mplog2(l)));
    1870         [ -  + ]:    4556861 :     if (!signe(x)) return NULL;
    1871                 :            :   }
    1872                 :    4685939 :   return x;
    1873                 :            : }
    1874                 :            : 
    1875                 :            : static GEN
    1876                 :    4685260 : mpexp_basecase(GEN x)
    1877                 :            : {
    1878                 :    4685260 :   pari_sp av = avma;
    1879                 :    4685260 :   long sh, l = realprec(x);
    1880                 :            :   GEN y, z;
    1881                 :            : 
    1882                 :    4685260 :   y = modlog2(x, &sh);
    1883         [ -  + ]:    4685260 :   if (!y) { avma = av; return real2n(sh, l); }
    1884                 :    4685260 :   z = addsr(1, exp1r_abs(y));
    1885         [ +  + ]:    4685260 :   if (signe(y) < 0) z = invr(z);
    1886         [ +  + ]:    4685260 :   if (sh) {
    1887                 :    4556182 :     shiftr_inplace(z, sh);
    1888         [ +  - ]:    4556182 :     if (realprec(z) > l) z = rtor(z, l); /* spurious precision increase */
    1889                 :            :   }
    1890                 :            : #ifdef DEBUG
    1891                 :            : {
    1892                 :            :   GEN t = mplog(z), u = divrr(subrr(x, t),x);
    1893                 :            :   if (signe(u) && expo(u) > 5-prec2nbits(minss(l,realprec(t))))
    1894                 :            :     pari_err_BUG("exp");
    1895                 :            : }
    1896                 :            : #endif
    1897                 :    4685260 :   return gerepileuptoleaf(av, z); /* NOT affrr, precision often increases */
    1898                 :            : }
    1899                 :            : 
    1900                 :            : GEN
    1901                 :    4716345 : mpexp(GEN x)
    1902                 :            : {
    1903                 :    4716345 :   const long s = 6; /*Initial steps using basecase*/
    1904                 :    4716345 :   long i, p, l = realprec(x), sh;
    1905                 :            :   GEN a, t, z;
    1906                 :            :   ulong mask;
    1907                 :            : 
    1908         [ +  + ]:    4716345 :   if (l <= maxss(EXPNEWTON_LIMIT, (1L<<s) + 2))
    1909                 :            :   {
    1910         [ +  + ]:    4715666 :     if (!signe(x)) return mpexp0(x);
    1911                 :    4684581 :     return mpexp_basecase(x);
    1912                 :            :   }
    1913                 :        679 :   z = cgetr(l); /* room for result */
    1914                 :        679 :   x = modlog2(x, &sh);
    1915         [ -  + ]:        679 :   if (!x) { avma = (pari_sp)(z+lg(z)); return real2n(sh, l); }
    1916                 :        679 :   constpi(l); /* precompute for later logr_abs() */
    1917                 :        679 :   mask = quadratic_prec_mask(prec2nbits(l)+BITS_IN_LONG);
    1918 [ +  + ][ +  + ]:       8826 :   for(i=0, p=1; i<s+TWOPOTBITS_IN_LONG; i++) { p <<= 1; if (mask & 1) p-=1; mask >>= 1; }
    1919                 :        679 :   a = mpexp_basecase(rtor(x, nbits2prec(p)));
    1920                 :        679 :   x = addrs(x,1);
    1921         [ -  + ]:        679 :   if (realprec(x) < l+EXTRAPRECWORD) x = rtor(x, l+EXTRAPRECWORD);
    1922                 :        679 :   a = rtor(a, l+EXTRAPRECWORD); /*append 0s */
    1923                 :        679 :   t = NULL;
    1924                 :            :   for(;;)
    1925                 :            :   {
    1926         [ -  + ]:        692 :     p <<= 1; if (mask & 1) p--;
    1927                 :        692 :     mask >>= 1;
    1928                 :        692 :     setprec(x, nbits2prec(p));
    1929                 :        692 :     setprec(a, nbits2prec(p));
    1930                 :        692 :     t = mulrr(a, subrr(x, logr_abs(a))); /* a (x - log(a)) */
    1931         [ +  + ]:        692 :     if (mask == 1) break;
    1932                 :         13 :     affrr(t, a); avma = (pari_sp)a;
    1933                 :         13 :   }
    1934                 :        679 :   affrr(t,z);
    1935         [ +  - ]:        679 :   if (sh) shiftr_inplace(z, sh);
    1936                 :    4716345 :   avma = (pari_sp)z; return z;
    1937                 :            : }
    1938                 :            : 
    1939                 :            : static long
    1940                 :       9975 : Qp_exp_prec(GEN x)
    1941                 :            : {
    1942                 :       9975 :   long k, e = valp(x), n = e + precp(x);
    1943                 :       9975 :   GEN p = gel(x,2);
    1944                 :       9975 :   int is2 = equaliu(p,2);
    1945 [ +  + ][ +  + ]:       9975 :   if (e < 1 || (e == 1 && is2)) return -1;
                 [ -  + ]
    1946         [ +  + ]:       9954 :   if (is2)
    1947                 :            :   {
    1948                 :       3017 :     n--; e--; k = n/e;
    1949         [ +  + ]:       3017 :     if (n%e == 0) k--;
    1950                 :            :   }
    1951                 :            :   else
    1952                 :            :   { /* e > 0, n > 0 */
    1953                 :       6937 :     GEN r, t = subis(p, 1);
    1954                 :       6937 :     k = itos(dvmdii(subis(muliu(t,n), 1), subis(muliu(t,e), 1), &r));
    1955         [ +  + ]:       6937 :     if (!signe(r)) k--;
    1956                 :            :   }
    1957                 :       9975 :   return k;
    1958                 :            : }
    1959                 :            : 
    1960                 :            : static GEN
    1961                 :      10801 : Qp_exp_safe(GEN x)
    1962                 :            : {
    1963                 :            :   long k;
    1964                 :            :   pari_sp av;
    1965                 :            :   GEN y;
    1966                 :            : 
    1967         [ +  + ]:      10801 :   if (gequal0(x)) return gaddgs(x,1);
    1968                 :       9919 :   k = Qp_exp_prec(x);
    1969         [ +  + ]:       9919 :   if (k < 0) return NULL;
    1970                 :       9912 :   av = avma;
    1971         [ +  + ]:      46837 :   for (y=gen_1; k; k--) y = gaddsg(1, gdivgs(gmul(y,x), k));
    1972                 :      10801 :   return gerepileupto(av, y);
    1973                 :            : }
    1974                 :            : 
    1975                 :            : GEN
    1976                 :      10339 : Qp_exp(GEN x)
    1977                 :            : {
    1978                 :      10339 :   GEN y = Qp_exp_safe(x);
    1979         [ +  + ]:      10339 :   if (!y) pari_err_DOMAIN("gexp(t_PADIC)","argument","",gen_0,x);
    1980                 :      10332 :   return y;
    1981                 :            : }
    1982                 :            : 
    1983                 :            : static GEN
    1984                 :         49 : cos_p(GEN x)
    1985                 :            : {
    1986                 :            :   long k;
    1987                 :            :   pari_sp av;
    1988                 :            :   GEN x2, y;
    1989                 :            : 
    1990         [ +  + ]:         49 :   if (gequal0(x)) return gaddgs(x,1);
    1991                 :         28 :   k = Qp_exp_prec(x);
    1992         [ +  + ]:         28 :   if (k < 0) return NULL;
    1993                 :         21 :   av = avma; x2 = gsqr(x);
    1994         [ -  + ]:         21 :   if (k & 1) k--;
    1995         [ +  + ]:        105 :   for (y=gen_1; k; k-=2)
    1996                 :            :   {
    1997                 :         84 :     GEN t = gdiv(gmul(y,x2), muluu(k, k-1));
    1998                 :         84 :     y = gsubsg(1, t);
    1999                 :            :   }
    2000                 :         49 :   return gerepileupto(av, y);
    2001                 :            : }
    2002                 :            : static GEN
    2003                 :         49 : sin_p(GEN x)
    2004                 :            : {
    2005                 :            :   long k;
    2006                 :            :   pari_sp av;
    2007                 :            :   GEN x2, y;
    2008                 :            : 
    2009         [ +  + ]:         49 :   if (gequal0(x)) return gcopy(x);
    2010                 :         28 :   k = Qp_exp_prec(x);
    2011         [ +  + ]:         28 :   if (k < 0) return NULL;
    2012                 :         21 :   av = avma; x2 = gsqr(x);
    2013         [ -  + ]:         21 :   if (k & 1) k--;
    2014         [ +  + ]:        105 :   for (y=gen_1; k; k-=2)
    2015                 :            :   {
    2016                 :         84 :     GEN t = gdiv(gmul(y,x2), muluu(k, k+1));
    2017                 :         84 :     y = gsubsg(1, t);
    2018                 :            :   }
    2019                 :         49 :   return gerepileupto(av, gmul(y, x));
    2020                 :            : }
    2021                 :            : 
    2022                 :            : static GEN
    2023                 :     713014 : cxexp(GEN x, long prec)
    2024                 :            : {
    2025                 :     713014 :   GEN r,p1,p2, y = cgetg(3,t_COMPLEX);
    2026                 :     713014 :   pari_sp av = avma, tetpil;
    2027                 :     713014 :   r = gexp(gel(x,1),prec);
    2028         [ -  + ]:     713014 :   if (gequal0(r)) { gel(y,1) = r; gel(y,2) = r; return y; }
    2029                 :     713014 :   gsincos(gel(x,2),&p2,&p1,prec);
    2030                 :     713014 :   tetpil = avma;
    2031                 :     713014 :   gel(y,1) = gmul(r,p1);
    2032                 :     713014 :   gel(y,2) = gmul(r,p2);
    2033                 :     713014 :   gerepilecoeffssp(av,tetpil,y+1,2);
    2034                 :     713014 :   return y;
    2035                 :            : }
    2036                 :            : 
    2037                 :            : /* given a t_SER x^v s(x), with s(0) != 0, return x^v(s - s(0)), shallow */
    2038                 :            : GEN
    2039                 :       2170 : serchop0(GEN s)
    2040                 :            : {
    2041                 :       2170 :   long i, l = lg(s);
    2042                 :            :   GEN y;
    2043         [ -  + ]:       2170 :   if (l == 2) return s;
    2044 [ +  + ][ -  + ]:       2170 :   if (l == 3 && isexactzero(gel(s,2))) return s;
    2045                 :       2170 :   y = cgetg(l, t_SER); y[1] = s[1];
    2046         [ +  + ]:      11683 :   gel(y,2) = gen_0; for (i=3; i <l; i++) gel(y,i) = gel(s,i);
    2047                 :       2170 :   return normalize(y);
    2048                 :            : }
    2049                 :            : 
    2050                 :            : static GEN
    2051                 :      19397 : serexp(GEN x, long prec)
    2052                 :            : {
    2053                 :            :   pari_sp av;
    2054                 :            :   long i,j,lx,ly,ex,mi;
    2055                 :            :   GEN p1,y,xd,yd;
    2056                 :            : 
    2057                 :      19397 :   ex = valp(x);
    2058         [ +  + ]:      19397 :   if (ex < 0) pari_err_DOMAIN("exp","valuation", "<", gen_0, x);
    2059         [ +  + ]:      19390 :   if (gequal0(x)) return gaddsg(1,x);
    2060                 :      18683 :   lx = lg(x);
    2061         [ +  + ]:      18683 :   if (ex)
    2062                 :            :   {
    2063                 :      17101 :     ly = lx+ex; y = cgetg(ly,t_SER);
    2064 [ +  + ][ +  + ]:      22015 :     mi = lx-1; while (mi>=3 && isrationalzero(gel(x,mi))) mi--;
    2065                 :      17101 :     mi += ex-2;
    2066                 :      17101 :     y[1] = evalsigne(1) | _evalvalp(0) | evalvarn(varn(x));
    2067                 :            :     /* zd[i] = coefficient of X^i in z */
    2068                 :      17101 :     xd = x+2-ex; yd = y+2; ly -= 2;
    2069                 :      17101 :     gel(yd,0) = gen_1;
    2070         [ +  + ]:      17311 :     for (i=1; i<ex; i++) gel(yd,i) = gen_0;
    2071         [ +  + ]:      57848 :     for (   ; i<ly; i++)
    2072                 :            :     {
    2073                 :      40747 :       av = avma; p1 = gen_0;
    2074         [ +  + ]:     123585 :       for (j=ex; j<=minss(i,mi); j++)
    2075                 :      82838 :         p1 = gadd(p1, gmulgs(gmul(gel(xd,j),gel(yd,i-j)),j));
    2076                 :      40747 :       gel(yd,i) = gerepileupto(av, gdivgs(p1,i));
    2077                 :            :     }
    2078                 :      17101 :     return y;
    2079                 :            :   }
    2080                 :       1582 :   av = avma;
    2081                 :      19390 :   return gerepileupto(av, gmul(gexp(gel(x,2),prec), serexp(serchop0(x),prec)));
    2082                 :            : }
    2083                 :            : 
    2084                 :            : GEN
    2085                 :    2289213 : gexp(GEN x, long prec)
    2086                 :            : {
    2087   [ +  +  +  + ]:    2289213 :   switch(typ(x))
    2088                 :            :   {
    2089                 :    1538139 :     case t_REAL: return mpexp(x);
    2090                 :     713014 :     case t_COMPLEX: return cxexp(x,prec);
    2091                 :         42 :     case t_PADIC: return Qp_exp(x);
    2092                 :            :     default:
    2093                 :            :     {
    2094                 :      38018 :       pari_sp av = avma;
    2095                 :            :       GEN y;
    2096         [ +  + ]:      38018 :       if (!(y = toser_i(x))) break;
    2097                 :      17815 :       return gerepileupto(av, serexp(y,prec));
    2098                 :            :     }
    2099                 :            :   }
    2100                 :    2289199 :   return trans_eval("exp",gexp,x,prec);
    2101                 :            : }
    2102                 :            : 
    2103                 :            : /********************************************************************/
    2104                 :            : /**                                                                **/
    2105                 :            : /**                           AGM(X, Y)                            **/
    2106                 :            : /**                                                                **/
    2107                 :            : /********************************************************************/
    2108                 :            : static int
    2109                 :    2695704 : agmr_gap(GEN a, GEN b, long L)
    2110                 :            : {
    2111                 :    2695704 :   GEN d = subrr(b, a);
    2112 [ +  + ][ +  + ]:    2695704 :   return (signe(d) && expo(d) - expo(b) >= L);
    2113                 :            : }
    2114                 :            : /* assume x > 0 */
    2115                 :            : static GEN
    2116                 :     195916 : agm1r_abs(GEN x)
    2117                 :            : {
    2118                 :     195916 :   long l = realprec(x), L = 5-prec2nbits(l);
    2119                 :     195916 :   GEN a1, b1, y = cgetr(l);
    2120                 :     195916 :   pari_sp av = avma;
    2121                 :            : 
    2122                 :     195916 :   a1 = addrr(real_1(l), x); shiftr_inplace(a1, -1);
    2123                 :     195916 :   b1 = sqrtr_abs(x);
    2124         [ +  + ]:    2695704 :   while (agmr_gap(a1,b1,L))
    2125                 :            :   {
    2126                 :    2499788 :     GEN a = a1;
    2127                 :    2499788 :     a1 = addrr(a,b1); shiftr_inplace(a1, -1);
    2128                 :    2499788 :     b1 = sqrtr_abs(mulrr(a,b1));
    2129                 :            :   }
    2130                 :     195916 :   affrr_fixlg(a1,y); avma = av; return y;
    2131                 :            : }
    2132                 :            : 
    2133                 :            : struct agmcx_gap_t { long L, ex, cnt; };
    2134                 :            : 
    2135                 :            : static void
    2136                 :      17594 : agmcx_init(GEN x, long *prec, struct agmcx_gap_t *S)
    2137                 :            : {
    2138                 :      17594 :   long l = precision(x);
    2139         [ +  + ]:      17594 :   if (l) *prec = l;
    2140                 :      17594 :   S->L = 1-prec2nbits(*prec);
    2141                 :      17594 :   S->cnt = 0;
    2142                 :      17594 :   S->ex = LONG_MAX;
    2143                 :      17594 : }
    2144                 :            : 
    2145                 :            : static long
    2146                 :      17594 : agmcx_a_b(GEN x, GEN *a1, GEN *b1, long prec)
    2147                 :            : {
    2148                 :      17594 :   long rotate = 0;
    2149         [ +  + ]:      17594 :   if (gsigne(real_i(x))<0)
    2150                 :            :   { /* Rotate by +/-Pi/2, so that the choice of the principal square
    2151                 :            :      * root gives the optimal AGM. So a1 = +/-I*a1, b1=sqrt(-x). */
    2152         [ +  + ]:        553 :     if (gsigne(imag_i(x))<0) { *a1=mulcxI(*a1);  rotate=-1; }
    2153                 :        301 :     else                     { *a1=mulcxmI(*a1); rotate=1; }
    2154                 :        553 :     x = gneg(x);
    2155                 :            :   }
    2156                 :      17594 :   *b1 = gsqrt(x, prec);
    2157                 :      17594 :   return rotate;
    2158                 :            : }
    2159                 :            : /* return 0 if we must stop the AGM loop (a=b or a ~ b), 1 otherwise */
    2160                 :            : static int
    2161                 :     302105 : agmcx_gap(GEN a, GEN b, struct agmcx_gap_t *S)
    2162                 :            : {
    2163                 :     302105 :   GEN d = gsub(b, a);
    2164                 :     302105 :   long ex = S->ex;
    2165                 :     302105 :   S->ex = gexpo(d);
    2166 [ +  + ][ -  + ]:     302105 :   if (gequal0(d) || S->ex - gexpo(b) < S->L) return 0;
    2167                 :            :   /* if (S->ex >= ex) we're no longer making progress; twice in a row */
    2168         [ +  + ]:     289720 :   if (S->ex < ex) S->cnt = 0;
    2169                 :            :   else
    2170         [ +  + ]:      10446 :     if (S->cnt++) return 0;
    2171                 :     302105 :   return 1;
    2172                 :            : }
    2173                 :            : static GEN
    2174                 :      17580 : agm1cx(GEN x, long prec)
    2175                 :            : {
    2176                 :            :   struct agmcx_gap_t S;
    2177                 :            :   GEN a1, b1;
    2178                 :      17580 :   pari_sp av = avma;
    2179                 :            :   long rotate;
    2180                 :      17580 :   agmcx_init(x, &prec, &S);
    2181                 :      17580 :   a1 = gtofp(gmul2n(gadd(real_1(prec), x), -1), prec);
    2182                 :      17580 :   rotate = agmcx_a_b(x, &a1, &b1, prec);
    2183         [ +  + ]:     302019 :   while (agmcx_gap(a1,b1,&S))
    2184                 :            :   {
    2185                 :     284439 :     GEN a = a1;
    2186                 :     284439 :     a1 = gmul2n(gadd(a,b1),-1);
    2187                 :     284439 :     b1 = gsqrt(gmul(a,b1), prec);
    2188                 :            :   }
    2189 [ +  + ][ +  + ]:      17580 :   if (rotate) a1 = rotate>0 ? mulcxI(a1):mulcxmI(a1);
    2190                 :      17580 :   return gerepilecopy(av,a1);
    2191                 :            : }
    2192                 :            : 
    2193                 :            : GEN
    2194                 :         14 : zellagmcx(GEN a0, GEN b0, GEN r, GEN t, long prec)
    2195                 :            : {
    2196                 :            :   struct agmcx_gap_t S;
    2197                 :         14 :   pari_sp av = avma;
    2198                 :         14 :   GEN x = gdiv(a0, b0), a1, b1;
    2199                 :            :   long rotate;
    2200                 :         14 :   agmcx_init(x, &prec, &S);
    2201                 :         14 :   a1 = gtofp(gmul2n(gadd(real_1(prec), x), -1), prec);
    2202                 :         14 :   r = gsqrt(gdiv(gmul(a1,gaddgs(r, 1)),gadd(r, x)), prec);
    2203                 :         14 :   t = gmul(r, t);
    2204                 :         14 :   rotate = agmcx_a_b(x, &a1, &b1, prec);
    2205         [ +  + ]:         86 :   while (agmcx_gap(a1,b1,&S))
    2206                 :            :   {
    2207                 :         72 :     GEN a = a1, b = b1;
    2208                 :         72 :     a1 = gmul2n(gadd(a,b),-1);
    2209                 :         72 :     b1 = gsqrt(gmul(a,b), prec);
    2210                 :         72 :     r = gsqrt(gdiv(gmul(a1,gaddgs(r, 1)),gadd(gmul(b, r), a )), prec);
    2211                 :         72 :     t = gmul(r, t);
    2212                 :            :   }
    2213 [ -  + ][ #  # ]:         14 :   if (rotate) a1 = rotate>0 ? mulcxI(a1):mulcxmI(a1);
    2214                 :         14 :   a1 = gmul(a1, b0);
    2215                 :         14 :   t = gatan(gdiv(a1,t), prec);
    2216                 :            :   /* send t to the fundamental domain if necessary */
    2217         [ +  + ]:         14 :   if (gsigne(real_i(t))<0) t = gadd(t, mppi(prec));
    2218                 :         14 :   return gerepileupto(av,gdiv(t,a1));
    2219                 :            : }
    2220                 :            : 
    2221                 :            : /* agm(1,x) */
    2222                 :            : static GEN
    2223                 :       1815 : agm1(GEN x, long prec)
    2224                 :            : {
    2225                 :            :   GEN p1, a1, b1, y;
    2226                 :            :   long l, l2, ep;
    2227                 :            :   pari_sp av;
    2228                 :            : 
    2229         [ -  + ]:       1815 :   if (gequal0(x)) return gcopy(x);
    2230   [ +  +  +  +  :       1815 :   switch(typ(x))
                      + ]
    2231                 :            :   {
    2232                 :            :     case t_INT:
    2233         [ +  + ]:         28 :       if (!is_pm1(x)) break;
    2234         [ +  + ]:         21 :       return (signe(x) > 0)? real_1(prec): real_0(prec);
    2235                 :            : 
    2236         [ +  + ]:        555 :     case t_REAL: return signe(x) > 0? agm1r_abs(x): agm1cx(x, prec);
    2237                 :            : 
    2238                 :            :     case t_COMPLEX:
    2239         [ -  + ]:       1106 :       if (gequal0(gel(x,2))) return agm1(gel(x,1), prec);
    2240                 :       1106 :       return agm1cx(x, prec);
    2241                 :            : 
    2242                 :            :     case t_PADIC:
    2243                 :         14 :       av = avma;
    2244                 :         14 :       a1 = x; b1 = gen_1; l = precp(x);
    2245                 :            :       do
    2246                 :            :       {
    2247                 :         28 :         GEN a = a1;
    2248                 :         28 :         a1 = gmul2n(gadd(a,b1),-1);
    2249                 :         28 :         a = gmul(a,b1);
    2250         [ +  + ]:         28 :         b1 = Qp_sqrt(a); if (!b1) pari_err_SQRTN("Qp_sqrt",a);
    2251                 :         21 :         p1 = gsub(b1,a1); ep = valp(p1)-valp(b1);
    2252         [ -  + ]:         21 :         if (ep<=0) { b1 = gneg_i(b1); p1 = gsub(b1,a1); ep=valp(p1)-valp(b1); }
    2253                 :            :       }
    2254 [ +  + ][ +  - ]:         21 :       while (ep<l && !gequal0(p1));
    2255                 :          7 :       return gerepilecopy(av,a1);
    2256                 :            : 
    2257                 :            :     default:
    2258         [ +  + ]:        112 :       av = avma; if (!(y = toser_i(x))) break;
    2259                 :          7 :       a1 = y; b1 = gen_1; l = lg(y)-2;
    2260                 :          7 :       l2 = 5-prec2nbits(prec);
    2261                 :            :       do
    2262                 :            :       {
    2263                 :         28 :         GEN a = a1;
    2264                 :         28 :         a1 = gmul2n(gadd(a,b1),-1);
    2265                 :         28 :         b1 = gsqrt(gmul(a,b1), prec);
    2266                 :         28 :         p1 = gsub(b1,a1); ep = valp(p1)-valp(b1);
    2267                 :            :       }
    2268         [ +  - ]:         21 :       while (ep<l && !gequal0(p1)
    2269 [ +  + ][ +  - ]:         49 :                   && (!isinexactreal(p1) || gexpo(p1) - gexpo(b1) >= l2));
                 [ #  # ]
    2270                 :          7 :       return gerepilecopy(av,a1);
    2271                 :            :   }
    2272                 :       1808 :   return trans_eval("agm",agm1,x,prec);
    2273                 :            : }
    2274                 :            : 
    2275                 :            : GEN
    2276                 :       1703 : agm(GEN x, GEN y, long prec)
    2277                 :            : {
    2278                 :            :   pari_sp av;
    2279         [ +  + ]:       1703 :   if (is_matvec_t(typ(y)))
    2280                 :            :   {
    2281         [ +  + ]:         14 :     if (is_matvec_t(typ(x))) pari_err_TYPE2("agm",x,y);
    2282                 :          7 :     swap(x, y);
    2283                 :            :   }
    2284         [ -  + ]:       1696 :   if (gequal0(y)) return gcopy(y);
    2285                 :       1696 :   av = avma;
    2286                 :       1696 :   return gerepileupto(av, gmul(y, agm1(gdiv(x,y), prec)));
    2287                 :            : }
    2288                 :            : 
    2289                 :            : /********************************************************************/
    2290                 :            : /**                                                                **/
    2291                 :            : /**                             LOG(X)                             **/
    2292                 :            : /**                                                                **/
    2293                 :            : /********************************************************************/
    2294                 :            : /* atanh(u/v) using binary splitting */
    2295                 :            : static GEN
    2296                 :       3746 : atanhQ_split(ulong u, ulong v, long prec)
    2297                 :            : {
    2298                 :            :   long i, nmax;
    2299                 :       3746 :   GEN u2 = sqru(u), v2 = sqru(v);
    2300                 :       3746 :   double d = ((double)v) / u;
    2301                 :            :   struct abpq_res R;
    2302                 :            :   struct abpq A;
    2303                 :            :   /* satisfies (2n+1) (v/u)^2n > 2^bitprec */
    2304                 :       3746 :   nmax = (long)(prec2nbits(prec) / (2*log2(d)));
    2305                 :       3746 :   abpq_init(&A, nmax);
    2306                 :       3746 :   A.a[0] = A.b[0] = gen_1;
    2307                 :       3746 :   A.p[0] = utoipos(u);
    2308                 :       3746 :   A.q[0] = utoipos(v);
    2309         [ +  + ]:     302992 :   for (i = 1; i <= nmax; i++)
    2310                 :            :   {
    2311                 :     299246 :     A.a[i] = gen_1;
    2312                 :     299246 :     A.b[i] = utoipos((i<<1)+1);
    2313                 :     299246 :     A.p[i] = u2;
    2314                 :     299246 :     A.q[i] = v2;
    2315                 :            :   }
    2316                 :       3746 :   abpq_sum(&R, 0, nmax, &A);
    2317                 :       3746 :   return rdivii(R.T, mulii(R.B,R.Q),prec);
    2318                 :            : }
    2319                 :            : /* log(2) = 10*atanh(1/17)+4*atanh(13/499) */
    2320                 :            : static GEN
    2321                 :       1873 : log2_split(long prec)
    2322                 :            : {
    2323                 :       1873 :   GEN u = atanhQ_split(1, 17, prec);
    2324                 :       1873 :   GEN v = atanhQ_split(13, 499, prec);
    2325                 :       1873 :   shiftr_inplace(v, 2);
    2326                 :       1873 :   return addrr(mulur(10, u), v);
    2327                 :            : }
    2328                 :            : #if 0 /* slower ! */
    2329                 :            : /* cf logagmr_abs(). Compute Pi/2agm(1, 4/2^n) ~ log(2^n) = n log(2) */
    2330                 :            : static GEN
    2331                 :            : log2_agm(long prec)
    2332                 :            : {
    2333                 :            :   long n = prec2nbits(prec) >> 1;
    2334                 :            :   GEN y = divrr(Pi2n(-1, prec), agm1r_abs( real2n(2 - n, prec) ));
    2335                 :            :   return divru(y, n);
    2336                 :            : }
    2337                 :            : #endif
    2338                 :            : GEN
    2339                 :    6379022 : constlog2(long prec)
    2340                 :            : {
    2341                 :            :   pari_sp av;
    2342                 :            :   GEN tmp;
    2343 [ +  + ][ +  + ]:    6379022 :   if (glog2 && realprec(glog2) >= prec) return glog2;
    2344                 :            : 
    2345                 :       1873 :   tmp = cgetr_block(prec);
    2346                 :       1873 :   av = avma;
    2347                 :       1873 :   affrr(log2_split(prec+EXTRAPRECWORD), tmp);
    2348                 :       1873 :   swap_clone(&glog2,tmp);
    2349                 :    6379022 :   avma = av; return glog2;
    2350                 :            : }
    2351                 :            : 
    2352                 :            : GEN
    2353                 :    6379022 : mplog2(long prec) { return rtor(constlog2(prec), prec); }
    2354                 :            : 
    2355                 :            : static GEN
    2356                 :     291421 : logagmr_abs(GEN q)
    2357                 :            : {
    2358                 :     291421 :   long prec = realprec(q), lim, e = expo(q);
    2359                 :            :   GEN z, y, Q, _4ovQ;
    2360                 :            :   pari_sp av;
    2361                 :            : 
    2362 [ +  + ][ +  + ]:     291421 :   if (absrnz_equal2n(q)) return e? mulsr(e, mplog2(prec)): real_0(prec);
    2363                 :     195403 :   z = cgetr(prec); av = avma; incrprec(prec);
    2364                 :     195403 :   lim = prec2nbits(prec) >> 1;
    2365                 :     195403 :   Q = rtor(q,prec);
    2366                 :     195403 :   shiftr_inplace(Q,lim-e); setsigne(Q,1);
    2367                 :            : 
    2368                 :     195403 :   _4ovQ = invr(Q); shiftr_inplace(_4ovQ, 2); /* 4/Q */
    2369                 :            :   /* Pi / 2agm(1, 4/Q) ~ log(Q), q = Q * 2^(e-lim) */
    2370                 :     195403 :   y = divrr(Pi2n(-1, prec), agm1r_abs(_4ovQ));
    2371                 :     195403 :   y = addrr(y, mulsr(e - lim, mplog2(prec)));
    2372                 :     291421 :   affrr_fixlg(y, z); avma = av; return z;
    2373                 :            : }
    2374                 :            : /*return log(|x|), assuming x != 0 */
    2375                 :            : GEN
    2376                 :    1830813 : logr_abs(GEN X)
    2377                 :            : {
    2378                 :            :   pari_sp ltop;
    2379                 :    1830813 :   long EX, L, m, k, a, b, l = realprec(X);
    2380                 :            :   GEN z, x, y;
    2381                 :            :   ulong u;
    2382                 :            :   double d;
    2383                 :            : 
    2384         [ +  + ]:    1830813 :   if (l > LOGAGM_LIMIT) return logagmr_abs(X);
    2385                 :            : 
    2386                 :            :  /* Assuming 1 < x < 2, we want delta = x-1, 1-x/2, 1-1/x, or 2/x-1 small.
    2387                 :            :   * We have 2/x-1 > 1-x/2, 1-1/x < x-1. So one should be choosing between
    2388                 :            :   * 1-1/x and 1-x/2 ( crossover sqrt(2), worse ~ 0.29 ). To avoid an inverse,
    2389                 :            :   * we choose between x-1 and 1-x/2 ( crossover 4/3, worse ~ 0.33 ) */
    2390                 :    1539392 :   EX = expo(X);
    2391                 :    1539392 :   u = uel(X,2);
    2392                 :    1539392 :   k = 2;
    2393         [ +  + ]:    1539392 :   if (u > (~0UL / 3) * 2) { /* choose 1-x/2 */
    2394                 :     851831 :     EX++; u = ~u;
    2395 [ +  + ][ +  + ]:     853998 :     while (!u && ++k < l) { u = uel(X,k); u = ~u; }
    2396                 :            :   } else { /* choose x - 1 */
    2397                 :     687561 :     u &= ~HIGHBIT; /* u - HIGHBIT, assuming HIGHBIT set */
    2398 [ +  + ][ +  + ]:    1224050 :     while (!u && ++k < l) u = uel(X,k);
    2399                 :            :   }
    2400 [ +  + ][ +  + ]:    1539392 :   if (k == l) return EX? mulsr(EX, mplog2(l)): real_0(l);
    2401         [ +  + ]:    1472543 :   z = cgetr(EX? l: l - (k-2)); ltop = avma;
    2402                 :            : 
    2403 [ +  + ][ +  + ]:    1472543 :   a = prec2nbits(k) + bfffo(u); /* ~ -log2 |1-x| */
         [ +  + ][ +  + ]
    2404                 :            :  /* Multiplication is quadratic in this range (l is small, otherwise we
    2405                 :            :   * use AGM). Set Y = x^(1/2^m), y = (Y - 1) / (Y + 1) and compute truncated
    2406                 :            :   * series sum y^(2k+1)/(2k+1): the costs is less than
    2407                 :            :   *    m b^2 + sum_{k <= n} ((2k+1) e + BITS_IN_LONG)^2
    2408                 :            :   * bit operations with |x-1| <  2^(1-a), |Y| < 2^(1-e) , m = e-a and b bits of
    2409                 :            :   * accuracy needed (+ BITS_IN_LONG since bit accuracies increase by
    2410                 :            :   * increments of BITS_IN_LONG), so
    2411                 :            :   * 4n^3/3 e^2 + n^2 2e BITS_IN_LONG+ n BITS_IN_LONG ~ m b^2, with n ~ b/2e
    2412                 :            :   * or b/6e + BITS_IN_LONG/2e + BITS_IN_LONG/2be ~ m
    2413                 :            :   *    B := (b / 6 + BITS_IN_LONG/2 + BITS_IN_LONG^2 / 2b) ~ m(m+a)
    2414                 :            :   *     m = min( -a/2 + sqrt(a^2/4 + B),  b - a )
    2415                 :            :   * NB: e ~ (b/6)^(1/2) as b -> oo */
    2416                 :    1472543 :   L = l+1;
    2417                 :    1472543 :   b = prec2nbits(L - (k-2)); /* take loss of accuracy into account */
    2418                 :            :   /* instead of the above pessimistic estimate for the cost of the sum, use
    2419                 :            :    * optimistic estimate (BITS_IN_LONG -> 0) */
    2420                 :    1472543 :   d = -a/2.; m = (long)(d + sqrt(d*d + b/6)); /* >= 0 */
    2421                 :            : 
    2422         [ +  + ]:    1472543 :   if (m > b-a) m = b-a;
    2423         [ +  + ]:    1472543 :   if (m < 0.2*a) m = 0; else L += nbits2extraprec(m);
    2424                 :    1472543 :   x = rtor(X,L);
    2425                 :    1472543 :   setsigne(x,1); shiftr_inplace(x,-EX);
    2426                 :            :   /* 2/3 < x < 4/3 */
    2427         [ +  + ]:    8761386 :   for (k=1; k<=m; k++) x = sqrtr_abs(x);
    2428                 :            : 
    2429                 :    1472543 :   y = divrr(subrs(x,1), addrs(x,1)); /* = (x-1) / (x+1), close to 0 */
    2430                 :    1472543 :   L = realprec(y); /* should be ~ l+1 - (k-2) */
    2431                 :            :   /* log(x) = log(1+y) - log(1-y) = 2 sum_{k odd} y^k / k
    2432                 :            :    * Truncate the sum at k = 2n+1, the remainder is
    2433                 :            :    *   2 sum_{k >= 2n+3} y^k / k < 2y^(2n+3) / (2n+3)(1-y) < y^(2n+3)
    2434                 :            :    * We want y^(2n+3) < y 2^(-prec2nbits(L)), hence
    2435                 :            :    *   n+1 > -prec2nbits(L) /-log_2(y^2) */
    2436                 :    1472543 :   d = -2*dbllog2r(y); /* ~ -log_2(y^2) */
    2437                 :    1472543 :   k = (long)(2*(prec2nbits(L) / d));
    2438                 :    1472543 :   k |= 1;
    2439         [ +  + ]:    1472543 :   if (k >= 3)
    2440                 :            :   {
    2441                 :    1471884 :     GEN S, T, y2 = sqrr(y), unr = real_1(L);
    2442                 :    1471884 :     pari_sp av = avma;
    2443                 :    1471884 :     long s = 0, incs = (long)d, l1 = nbits2prec((long)d);
    2444                 :    1471884 :     S = x;
    2445                 :    1471884 :     setprec(S,  l1);
    2446                 :    1471884 :     setprec(unr,l1); affrr(divru(unr,k), S); /* destroy x, not needed anymore */
    2447                 :    1471884 :     for (k -= 2;; k -= 2) /* k = 2n+1, ..., 1 */
    2448                 :            :     { /* S = y^(2n+1-k)/(2n+1) + ... + 1 / k */
    2449                 :   26582888 :       setprec(y2, l1); T = mulrr(S,y2);
    2450         [ +  + ]:   26582888 :       if (k == 1) break;
    2451                 :            : 
    2452         [ -  + ]:   25111004 :       l1 += dvmdsBIL(s + incs, &s); if (l1>L) l1=L;
    2453                 :   25111004 :       setprec(S, l1);
    2454                 :   25111004 :       setprec(unr,l1);
    2455                 :   25111004 :       affrr(addrr(divru(unr, k), T), S); avma = av;
    2456                 :   25111004 :     }
    2457                 :            :     /* k = 1 special-cased for eficiency */
    2458                 :    1471884 :     y = mulrr(y, addsr(1,T)); /* = log(X)/2 */
    2459                 :            :   }
    2460                 :    1472543 :   shiftr_inplace(y, m + 1);
    2461         [ +  + ]:    1472543 :   if (EX) y = addrr(y, mulsr(EX, mplog2(l+1)));
    2462                 :    1830813 :   affrr_fixlg(y, z); avma = ltop; return z;
    2463                 :            : }
    2464                 :            : 
    2465                 :            : /* assume Im(q) != 0 and precision(q) >= prec. Compute log(q) with accuracy
    2466                 :            :  * prec [disregard input accuracy] */
    2467                 :            : GEN
    2468                 :      16432 : logagmcx(GEN q, long prec)
    2469                 :            : {
    2470                 :      16432 :   GEN z = cgetc(prec), y, Q, a, b;
    2471                 :            :   long lim, e, ea, eb;
    2472                 :      16432 :   pari_sp av = avma;
    2473                 :      16432 :   int neg = 0;
    2474                 :            : 
    2475                 :      16432 :   incrprec(prec);
    2476         [ +  + ]:      16432 :   if (gsigne(gel(q,1)) < 0) { q = gneg(q); neg = 1; }
    2477                 :      16432 :   lim = prec2nbits(prec) >> 1;
    2478                 :      16432 :   Q = gtofp(q, prec);
    2479                 :      16432 :   a = gel(Q,1);
    2480                 :      16432 :   b = gel(Q,2);
    2481         [ -  + ]:      16432 :   if (gequal0(a)) {
    2482                 :          0 :     affrr_fixlg(logr_abs(b), gel(z,1));
    2483                 :          0 :     y = Pi2n(-1, prec);
    2484         [ #  # ]:          0 :     if (signe(b) < 0) setsigne(y, -1);
    2485                 :          0 :     affrr_fixlg(y, gel(z,2)); avma = av; return z;
    2486                 :            :   }
    2487                 :      16432 :   ea = expo(a);
    2488                 :      16432 :   eb = expo(b);
    2489         [ +  + ]:      16432 :   e = ea <= eb ? lim - eb : lim - ea;
    2490                 :      16432 :   shiftr_inplace(a, e);
    2491                 :      16432 :   shiftr_inplace(b, e);
    2492                 :            : 
    2493                 :            :   /* Pi / 2agm(1, 4/Q) ~ log(Q), q = Q * 2^e */
    2494                 :      16432 :   y = gdiv(Pi2n(-1, prec), agm1cx( gdivsg(4, Q), prec ));
    2495                 :      16432 :   a = gel(y,1);
    2496                 :      16432 :   b = gel(y,2);
    2497                 :      16432 :   a = addrr(a, mulsr(-e, mplog2(prec)));
    2498         [ +  + ]:      16432 :   if (realprec(a) <= LOWDEFAULTPREC) a = real_0_bit(expo(a));
    2499         [ +  + ]:      17438 :   if (neg) b = gsigne(b) <= 0? gadd(b, mppi(prec))
    2500         [ +  + ]:       1006 :                              : gsub(b, mppi(prec));
    2501                 :      16432 :   affrr_fixlg(a, gel(z,1));
    2502                 :      16432 :   affrr_fixlg(b, gel(z,2)); avma = av; return z;
    2503                 :            : }
    2504                 :            : 
    2505                 :            : GEN
    2506                 :     245903 : mplog(GEN x)
    2507                 :            : {
    2508         [ -  + ]:     245903 :   if (signe(x)<=0) pari_err_DOMAIN("mplog", "argument", "<=", gen_0, x);
    2509                 :     245903 :   return logr_abs(x);
    2510                 :            : }
    2511                 :            : 
    2512                 :            : /* pe = p^e, p prime, 0 < x < pe a t_INT coprime to p. Return the (p-1)-th
    2513                 :            :  * root of 1 in (Z/pe)^* congruent to x mod p, resp x mod 4 if p = 2.
    2514                 :            :  * Simplified form of Zp_sqrtnlift: 1/(p-1) is trivial to compute */
    2515                 :            : GEN
    2516                 :          7 : Zp_teichmuller(GEN x, GEN p, long e, GEN pe)
    2517                 :            : {
    2518                 :            :   GEN q, z, p1;
    2519                 :            :   pari_sp av;
    2520                 :            :   ulong mask;
    2521 [ -  + ][ #  # ]:          7 :   if (equaliu(p,2)) return (mod4(x) & 2)? subiu(pe,1): gen_1;
    2522         [ -  + ]:          7 :   if (e == 1) return icopy(x);
    2523                 :          7 :   av = avma;
    2524                 :          7 :   p1 = subiu(p, 1);
    2525                 :          7 :   mask = quadratic_prec_mask(e);
    2526                 :          7 :   q = p; z = remii(x, p);
    2527         [ +  + ]:         35 :   while (mask > 1)
    2528                 :            :   { /* Newton iteration solving z^{1 - p} = 1, z = x (mod p) */
    2529                 :         28 :     GEN w, t, qold = q;
    2530         [ +  + ]:         28 :     if (mask <= 3) /* last iteration */
    2531                 :          7 :       q = pe;
    2532                 :            :     else
    2533                 :            :     {
    2534                 :         21 :       q = sqri(q);
    2535         [ +  + ]:         21 :       if (mask & 1) q = diviiexact(q, p);
    2536                 :            :     }
    2537                 :         28 :     mask >>= 1;
    2538                 :            :     /* q <= qold^2 */
    2539         [ +  + ]:         28 :     if (lgefint(q) == 3)
    2540                 :            :     {
    2541                 :         20 :       ulong Z = uel(z,2), Q = uel(q,2), P1 = uel(p1,2);
    2542                 :         20 :       ulong W = (Q-1) / P1; /* -1/(p-1) + O(qold) */
    2543                 :         20 :       ulong T = Fl_mul(W, Fl_powu(Z,P1,Q) - 1, Q);
    2544                 :         20 :       Z = Fl_mul(Z, 1 + T, Q);
    2545                 :         20 :       z = utoi(Z);
    2546                 :            :     }
    2547                 :            :     else
    2548                 :            :     {
    2549                 :          8 :       w = diviiexact(addsi(-1,qold),p1); /* -1/(p-1) + O(qold) */
    2550                 :          8 :       t = Fp_mul(w, subis(Fp_pow(z,p1,q), 1), q);
    2551                 :          8 :       z = Fp_mul(z, addsi(1,t), q);
    2552                 :            :     }
    2553                 :            :   }
    2554                 :          7 :   return gerepileuptoint(av, z);
    2555                 :            : }
    2556                 :            : 
    2557                 :            : GEN
    2558                 :          7 : teich(GEN x)
    2559                 :            : {
    2560                 :            :   GEN p, q, y, z;
    2561                 :            :   long n;
    2562                 :            : 
    2563         [ -  + ]:          7 :   if (typ(x)!=t_PADIC) pari_err_TYPE("teichmuller",x);
    2564                 :          7 :   z = gel(x,4);
    2565         [ -  + ]:          7 :   if (!signe(z)) return gcopy(x);
    2566                 :          7 :   p = gel(x,2);
    2567                 :          7 :   q = gel(x,3);
    2568                 :          7 :   n = precp(x);
    2569                 :          7 :   y = cgetg(5,t_PADIC);
    2570                 :          7 :   y[1] = evalprecp(n) | _evalvalp(0);
    2571                 :          7 :   gel(y,2) = icopy(p);
    2572                 :          7 :   gel(y,3) = icopy(q);
    2573                 :          7 :   gel(y,4) = Zp_teichmuller(z, p, n, q);
    2574                 :          7 :   return y;
    2575                 :            : }
    2576                 :            : 
    2577                 :            : GEN
    2578                 :    1681027 : glog(GEN x, long prec)
    2579                 :            : {
    2580                 :            :   pari_sp av, tetpil;
    2581                 :            :   GEN y, p1;
    2582                 :            :   long l;
    2583                 :            : 
    2584   [ +  +  +  +  :    1681027 :   switch(typ(x))
                      + ]
    2585                 :            :   {
    2586                 :            :     case t_REAL:
    2587         [ +  + ]:    1071856 :       if (signe(x) >= 0)
    2588                 :            :       {
    2589         [ +  + ]:     836251 :         if (!signe(x)) pari_err_DOMAIN("log", "argument", "=", gen_0, x);
    2590                 :     836237 :         return logr_abs(x);
    2591                 :            :       }
    2592                 :     235605 :       retmkcomplex(logr_abs(x), mppi(realprec(x)));
    2593                 :            : 
    2594                 :            :     case t_FRAC:
    2595                 :            :     {
    2596                 :            :       GEN a, b;
    2597                 :            :       long e1, e2;
    2598                 :      93090 :       av = avma;
    2599                 :      93090 :       a = gel(x,1);
    2600                 :      93090 :       b = gel(x,2);
    2601                 :      93090 :       e1 = expi(subii(a,b)); e2 = expi(b);
    2602         [ +  + ]:      93090 :       if (e2 > e1) prec += nbits2nlong(e2 - e1);
    2603                 :      93090 :       x = fractor(x, prec);
    2604                 :      93090 :       return gerepileupto(av, glog(x, prec));
    2605                 :            :     }
    2606                 :            :     case t_COMPLEX:
    2607         [ +  + ]:     224039 :       if (ismpzero(gel(x,2))) return glog(gel(x,1), prec);
    2608         [ +  + ]:     221626 :       if (ismpzero(gel(x,1)))
    2609                 :            :       {
    2610                 :       2667 :         GEN a = gel(x,2), b;
    2611                 :       2667 :         av = avma; b = Pi2n(-1,prec);
    2612         [ +  + ]:       2667 :         if (gsigne(a) < 0) { setsigne(b, -1); a = gabs(a,prec); }
    2613         [ +  + ]:       2667 :         a = isint1(a) ? gen_0: glog(a,prec);
    2614                 :       2667 :         return gerepilecopy(av, mkcomplex(a, b));
    2615                 :            :       }
    2616         [ +  + ]:     218959 :       l = precision(x); if (l > prec) prec = l;
    2617         [ +  + ]:     218959 :       if (prec >= LOGAGMCX_LIMIT) return logagmcx(x, prec);
    2618                 :     202723 :       y = cgetg(3,t_COMPLEX);
    2619                 :     202723 :       gel(y,2) = garg(x,prec);
    2620                 :     202723 :       av = avma; p1 = glog(cxnorm(x),prec); tetpil = avma;
    2621                 :     202723 :       gel(y,1) = gerepile(av,tetpil,gmul2n(p1,-1)); return y;
    2622                 :            : 
    2623                 :        231 :     case t_PADIC: return Qp_log(x);
    2624                 :            :     default:
    2625         [ +  + ]:     291811 :       av = avma; if (!(y = toser_i(x))) break;
    2626         [ -  + ]:         21 :       if (!signe(y)) pari_err_DOMAIN("log", "argument", "=", gen_0, x);
    2627         [ +  + ]:         21 :       if (valp(y)) pari_err_DOMAIN("log", "series valuation", "!=", gen_0, x);
    2628                 :         14 :       p1 = integser(gdiv(derivser(y), y)); /* log(y)' = y'/y */
    2629         [ +  + ]:         14 :       if (!gequal1(gel(y,2))) p1 = gadd(p1, glog(gel(y,2),prec));
    2630                 :         14 :       return gerepileupto(av, p1);
    2631                 :            :   }
    2632                 :    1680992 :   return trans_eval("log",glog,x,prec);
    2633                 :            : }
    2634                 :            : /********************************************************************/
    2635                 :            : /**                                                                **/
    2636                 :            : /**                        SINE, COSINE                            **/
    2637                 :            : /**                                                                **/
    2638                 :            : /********************************************************************/
    2639                 :            : 
    2640                 :            : /* Reduce x0 mod Pi/2 to x in [-Pi/4, Pi/4]. Return cos(x)-1 */
    2641                 :            : static GEN
    2642                 :    4303492 : mpcosm1(GEN x, long *ptmod8)
    2643                 :            : {
    2644                 :    4303492 :   long a = expo(x), l = realprec(x), b, L, i, n, m, B;
    2645                 :            :   GEN y, p2, x2;
    2646                 :            :   double d;
    2647                 :            : 
    2648                 :    4303492 :   n = 0;
    2649         [ +  + ]:    4303492 :   if (a >= 0)
    2650                 :            :   {
    2651                 :            :     long p;
    2652                 :            :     GEN q;
    2653         [ +  + ]:    3103814 :     if (a > 30)
    2654                 :            :     {
    2655                 :       4464 :       GEN z, pitemp = Pi2n(-2, nbits2prec(a + 32));
    2656                 :       4464 :       z = addrr(x,pitemp); /* = x + Pi/4 */
    2657         [ -  + ]:       4464 :       if (expo(z) >= bit_prec(z) + 3) pari_err_PREC("mpcosm1");
    2658                 :       4464 :       shiftr_inplace(pitemp, 1);
    2659                 :       4464 :       q = floorr( divrr(z,pitemp) ); /* round ( x / (Pi/2) ) */
    2660                 :       4464 :       p = l+EXTRAPRECWORD; x = rtor(x,p);
    2661                 :            :     } else {
    2662                 :    3099350 :       q = stoi((long)floor(rtodbl(x) / (PI/2) + 0.5));
    2663                 :    3099350 :       p = l;
    2664                 :            :     }
    2665         [ +  - ]:    3103814 :     if (signe(q))
    2666                 :            :     {
    2667                 :    3103814 :       x = subrr(x, mulir(q, Pi2n(-1,p))); /* x mod Pi/2  */
    2668                 :    3103814 :       a = expo(x);
    2669 [ +  + ][ -  + ]:    3103814 :       if (!signe(x) && a >= 0) pari_err_PREC("mpcosm1");
    2670 [ +  + ][ +  + ]:    3103814 :       n = mod4(q); if (n && signe(q) < 0) n = 4 - n;
    2671                 :            :     }
    2672                 :            :   }
    2673                 :            :   /* a < 0 */
    2674         [ +  + ]:    4303492 :   b = signe(x); *ptmod8 = (b < 0)? 4 + n: n;
    2675         [ +  + ]:    4303492 :   if (!b) return real_0_bit((expo(x)<<1) - 1);
    2676                 :            : 
    2677                 :    4133022 :   b = prec2nbits(l);
    2678         [ +  + ]:    4133022 :   if (b + (a<<1) <= 0) {
    2679                 :     206224 :     y = sqrr(x); shiftr_inplace(y, -1); setsigne(y, -1);
    2680                 :     206224 :     return y;
    2681                 :            :   }
    2682                 :            : 
    2683                 :    3926798 :   y = cgetr(l);
    2684                 :    3926798 :   B = b/6 + BITS_IN_LONG + (BITS_IN_LONG*BITS_IN_LONG/2)/ b;
    2685                 :    3926798 :   d = a/2.; m = (long)(d + sqrt(d*d + B)); /* >= 0 ,*/
    2686         [ +  + ]:    3926798 :   if (m < (-a) * 0.1) m = 0; /* not worth it */
    2687                 :    3926798 :   L = l + nbits2extraprec(m);
    2688                 :            : 
    2689                 :    3926798 :   b += m;
    2690                 :    3926798 :   d = 2.0 * (m-dbllog2r(x)-1/LOG2); /* ~ 2( - log_2 Y - 1/log(2) ) */
    2691                 :    3926798 :   n = (long)(b / d);
    2692         [ +  + ]:    3926798 :   if (n > 1)
    2693                 :    3865882 :     n = (long)(b / (d + log2((double)n+1))); /* log~constant in small ranges */
    2694         [ +  + ]:    7980748 :   while (n*(d+log2((double)n+1)) < b) n++; /* expect few corrections */
    2695                 :            : 
    2696                 :            :  /* Multiplication is quadratic in this range (l is small, otherwise we
    2697                 :            :   * use logAGM + Newton). Set Y = 2^(-e-a) x, compute truncated series
    2698                 :            :   * sum Y^2k/(2k)!: this costs roughly
    2699                 :            :   *   m b^2 + sum_{k <= n} (2k e + BITS_IN_LONG)^2
    2700                 :            :   *   ~ floor(b/2e) b^2 / 3  + m b^2
    2701                 :            :   * bit operations with |x| <  2^(1+a), |Y| < 2^(1-e) , m = e+a and b bits of
    2702                 :            :   * accuracy needed, so
    2703                 :            :   *    B := ( b / 6 + BITS_IN_LONG + BITS_IN_LONG^2 / 2b) ~ m(m-a)
    2704                 :            :   * we want b ~ 6 m (m-a) or m~b+a hence
    2705                 :            :   *     m = min( a/2 + sqrt(a^2/4 + b/6),  b/2 + a )
    2706                 :            :   * NB1: e ~ (b/6)^(1/2) or b/2.
    2707                 :            :   * NB2: We use b/4 instead of b/6 in the formula above: hand-optimized...
    2708                 :            :   *
    2709                 :            :   * Truncate the sum at k = n (>= 1), the remainder is
    2710                 :            :   * < sum_{k >= n+1} Y^2k / 2k! < Y^(2n+2) / (2n+2)!(1-Y^2) < Y^(2n+2)/(2n+1)!
    2711                 :            :   * We want ... <= Y^2 2^-b, hence -2n log_2 |Y| + log_2 (2n+1)! >= b
    2712                 :            :   *   log n! ~ (n + 1/2) log(n+1) - (n+1) + log(2Pi)/2,
    2713                 :            :   * error bounded by 1/6(n+1) <= 1/12. Finally, we want
    2714                 :            :   * 2n (-1/log(2) - log_2 |Y| + log_2(2n+2)) >= b  */
    2715                 :    3926798 :   x = rtor(x, L); shiftr_inplace(x, -m); setsigne(x, 1);
    2716                 :    3926798 :   x2 = sqrr(x);
    2717         [ -  + ]:    3926798 :   if (n == 1) { p2 = x2; shiftr_inplace(p2, -1); setsigne(p2, -1); } /*-Y^2/2*/
    2718                 :            :   else
    2719                 :            :   {
    2720                 :    3926798 :     GEN unr = real_1(L);
    2721                 :            :     pari_sp av;
    2722                 :    3926798 :     long s = 0, l1 = nbits2prec((long)(d + n + 16));
    2723                 :            : 
    2724                 :    3926798 :     p2 = cgetr(L); av = avma;
    2725         [ +  + ]:   25700398 :     for (i=n; i>=2; i--)
    2726                 :            :     {
    2727                 :            :       GEN p1;
    2728                 :   21773600 :       setprec(x2,l1); p1 = divrunu(x2, 2*i-1);
    2729         [ +  + ]:   21773600 :       l1 += dvmdsBIL(s - expo(p1), &s); if (l1>L) l1=L;
    2730         [ +  + ]:   21773600 :       if (i != n) p1 = mulrr(p1,p2);
    2731                 :   21773600 :       setprec(unr,l1); p1 = addrr_sign(unr,1, p1,-signe(p1));
    2732                 :   21773600 :       setprec(p2,l1); affrr(p1,p2); avma = av;
    2733                 :            :     }
    2734                 :    3926798 :     shiftr_inplace(p2, -1); togglesign(p2); /* p2 := -p2/2 */
    2735                 :    3926798 :     setprec(x2,L); p2 = mulrr(x2,p2);
    2736                 :            :   }
    2737                 :            :   /* Now p2 = sum {1<= i <=n} (-1)^i x^(2i) / (2i)! ~ cos(x) - 1 */
    2738         [ +  + ]:   35649468 :   for (i=1; i<=m; i++)
    2739                 :            :   { /* p2 = cos(x)-1 --> cos(2x)-1 */
    2740                 :   31722670 :     p2 = mulrr(p2, addsr(2,p2));
    2741                 :   31722670 :     shiftr_inplace(p2, 1);
    2742         [ +  + ]:   31722670 :     if ((i & 31) == 0) p2 = gerepileuptoleaf((pari_sp)y, p2);
    2743                 :            :   }
    2744                 :    4303492 :   affrr_fixlg(p2,y); return y;
    2745                 :            : }
    2746                 :            : 
    2747                 :            : /* sqrt (|1 - (1+x)^2|) = sqrt(|x*(x+2)|). Sends cos(x)-1 to |sin(x)| */
    2748                 :            : static GEN
    2749                 :    2734874 : mpaut(GEN x)
    2750                 :            : {
    2751                 :    2734874 :   pari_sp av = avma;
    2752                 :    2734874 :   GEN t = mulrr(x, addsr(2,x)); /* != 0 */
    2753         [ +  + ]:    2734874 :   if (!signe(t)) return real_0_bit(expo(t) >> 1);
    2754                 :    2734874 :   return gerepileuptoleaf(av, sqrtr_abs(t));
    2755                 :            : }
    2756                 :            : 
    2757                 :            : /********************************************************************/
    2758                 :            : /**                            COSINE                              **/
    2759                 :            : /********************************************************************/
    2760                 :            : 
    2761                 :            : GEN
    2762                 :    2758364 : mpcos(GEN x)
    2763                 :            : {
    2764                 :            :   long mod8;
    2765                 :            :   pari_sp av;
    2766                 :            :   GEN y,p1;
    2767                 :            : 
    2768         [ +  + ]:    2758364 :   if (!signe(x)) {
    2769                 :         32 :     long l = nbits2prec(-expo(x));
    2770         [ -  + ]:         32 :     if (l < LOWDEFAULTPREC) l = LOWDEFAULTPREC;
    2771                 :         32 :     return real_1(l);
    2772                 :            :   }
    2773                 :            : 
    2774                 :    2758332 :   av = avma; p1 = mpcosm1(x,&mod8);
    2775   [ +  +  +  + ]:    2758332 :   switch(mod8)
    2776                 :            :   {
    2777                 :     770468 :     case 0: case 4: y = addsr(1,p1); break;
    2778                 :     709785 :     case 1: case 7: y = mpaut(p1); togglesign(y); break;
    2779                 :     657208 :     case 2: case 6: y = subsr(-1,p1); break;
    2780                 :     620871 :     default:        y = mpaut(p1); break; /* case 3: case 5: */
    2781                 :            :   }
    2782                 :    2758364 :   return gerepileuptoleaf(av, y);
    2783                 :            : }
    2784                 :            : 
    2785                 :            : /* convert INT or FRAC to REAL, which is later reduced mod 2Pi : avoid
    2786                 :            :  * cancellation */
    2787                 :            : static GEN
    2788                 :       2569 : tofp_safe(GEN x, long prec)
    2789                 :            : {
    2790         [ +  + ]:       5138 :   return (typ(x) == t_INT || gexpo(x) > 0)? gadd(x, real_0(prec))
    2791         [ +  + ]:       5138 :                                           : fractor(x, prec);
    2792                 :            : }
    2793                 :            : 
    2794                 :            : GEN
    2795                 :     168172 : gcos(GEN x, long prec)
    2796                 :            : {
    2797                 :            :   pari_sp av;
    2798                 :            :   GEN r, u, v, y, u1, v1;
    2799                 :            :   long i;
    2800                 :            : 
    2801   [ +  +  +  +  :     168172 :   switch(typ(x))
                      + ]
    2802                 :            :   {
    2803                 :     167325 :     case t_REAL: return mpcos(x);
    2804                 :            :     case t_COMPLEX:
    2805         [ +  + ]:         63 :       if (isintzero(gel(x,1))) return gcosh(gel(x,2), prec);
    2806         [ -  + ]:         56 :       i = precision(x); if (!i) i = prec;
    2807                 :         56 :       y = cgetc(i); av = avma;
    2808                 :         56 :       r = gexp(gel(x,2),prec);
    2809                 :         56 :       v1 = gmul2n(addrr(invr(r),r), -1); /* = cos(I*Im(x)) */
    2810                 :         56 :       u1 = subrr(v1, r); /* = - I*sin(I*Im(x)) */
    2811                 :         56 :       gsincos(gel(x,1),&u,&v,prec);
    2812                 :         56 :       affrr_fixlg(gmul(v1,v), gel(y,1));
    2813                 :         56 :       affrr_fixlg(gmul(u1,u), gel(y,2)); avma = av; return y;
    2814                 :            : 
    2815                 :            :     case t_INT: case t_FRAC:
    2816                 :        707 :       y = cgetr(prec); av = avma;
    2817                 :        707 :       affrr_fixlg(mpcos(tofp_safe(x,prec)), y); avma = av; return y;
    2818                 :            : 
    2819                 :         49 :     case t_PADIC: y = cos_p(x);
    2820         [ +  + ]:         49 :       if (!y) pari_err_DOMAIN("gcos(t_PADIC)","argument","",gen_0,x);
    2821                 :         42 :       return y;
    2822                 :            : 
    2823                 :            :     default:
    2824         [ +  + ]:         28 :       av = avma; if (!(y = toser_i(x))) break;
    2825         [ -  + ]:         21 :       if (gequal0(y)) return gerepileupto(av, gaddsg(1,y));
    2826         [ +  + ]:         21 :       if (valp(y) < 0)
    2827                 :          7 :         pari_err_DOMAIN("cos","valuation", "<", gen_0, x);
    2828                 :         14 :       gsincos(y,&u,&v,prec);
    2829                 :         14 :       return gerepilecopy(av,v);
    2830                 :            :   }
    2831                 :     168158 :   return trans_eval("cos",gcos,x,prec);
    2832                 :            : }
    2833                 :            : /********************************************************************/
    2834                 :            : /**                             SINE                               **/
    2835                 :            : /********************************************************************/
    2836                 :            : 
    2837                 :            : GEN
    2838                 :     741345 : mpsin(GEN x)
    2839                 :            : {
    2840                 :            :   long mod8;
    2841                 :            :   pari_sp av;
    2842                 :            :   GEN y,p1;
    2843                 :            : 
    2844         [ +  + ]:     741345 :   if (!signe(x)) return real_0_bit(expo(x));
    2845                 :            : 
    2846                 :     740914 :   av = avma; p1 = mpcosm1(x,&mod8);
    2847   [ +  +  +  + ]:     740914 :   switch(mod8)
    2848                 :            :   {
    2849                 :     491965 :     case 0: case 6: y=mpaut(p1); break;
    2850                 :      71245 :     case 1: case 5: y=addsr(1,p1); break;
    2851                 :     108007 :     case 2: case 4: y=mpaut(p1); togglesign(y); break;
    2852                 :      69697 :     default:        y=subsr(-1,p1); break; /* case 3: case 7: */
    2853                 :            :   }
    2854                 :     741345 :   return gerepileuptoleaf(av, y);
    2855                 :            : }
    2856                 :            : 
    2857                 :            : GEN
    2858                 :     741660 : gsin(GEN x, long prec)
    2859                 :            : {
    2860                 :            :   pari_sp av;
    2861                 :            :   GEN r, u, v, y, v1, u1;
    2862                 :            :   long i;
    2863                 :            : 
    2864   [ +  +  +  +  :     741660 :   switch(typ(x))
                      + ]
    2865                 :            :   {
    2866                 :     739665 :     case t_REAL: return mpsin(x);
    2867                 :            :     case t_COMPLEX:
    2868         [ +  + ]:        182 :       if (isintzero(gel(x,1))) retmkcomplex(gen_0,gsinh(gel(x,2),prec));
    2869         [ -  + ]:        175 :       i = precision(x); if (!i) i = prec;
    2870                 :        175 :       y = cgetc(i); av = avma;
    2871                 :        175 :       r = gexp(gel(x,2),prec);
    2872                 :        175 :       v1 = gmul2n(addrr(invr(r),r), -1); /* = cos(I*Im(x)) */
    2873                 :        175 :       u1 = subrr(r, v1); /* = I*sin(I*Im(x)) */
    2874                 :        175 :       gsincos(gel(x,1),&u,&v,prec);
    2875                 :        175 :       affrr_fixlg(gmul(v1,u), gel(y,1));
    2876                 :        175 :       affrr_fixlg(gmul(u1,v), gel(y,2)); avma = av; return y;
    2877                 :            : 
    2878                 :            :     case t_INT: case t_FRAC:
    2879                 :       1680 :       y = cgetr(prec); av = avma;
    2880                 :       1680 :       affrr_fixlg(mpsin(tofp_safe(x,prec)), y); avma = av; return y;
    2881                 :            : 
    2882                 :         49 :     case t_PADIC: y = sin_p(x);
    2883         [ +  + ]:         49 :       if (!y) pari_err_DOMAIN("gsin(t_PADIC)","argument","",gen_0,x);
    2884                 :         42 :       return y;
    2885                 :            : 
    2886                 :            :     default:
    2887         [ +  + ]:         84 :       av = avma; if (!(y = toser_i(x))) break;
    2888         [ -  + ]:         77 :       if (gequal0(y)) return gerepilecopy(av, y);
    2889         [ +  + ]:         77 :       if (valp(y) < 0)
    2890                 :          7 :         pari_err_DOMAIN("sin","valuation", "<", gen_0, x);
    2891                 :         70 :       gsincos(y,&u,&v,prec);
    2892                 :         70 :       return gerepilecopy(av,u);
    2893                 :            :   }
    2894                 :     741646 :   return trans_eval("sin",gsin,x,prec);
    2895                 :            : }
    2896                 :            : /********************************************************************/
    2897                 :            : /**                       SINE, COSINE together                    **/
    2898                 :            : /********************************************************************/
    2899                 :            : 
    2900                 :            : void
    2901                 :     828299 : mpsincos(GEN x, GEN *s, GEN *c)
    2902                 :            : {
    2903                 :            :   long mod8;
    2904                 :            :   pari_sp av, tetpil;
    2905                 :            :   GEN p1, *gptr[2];
    2906                 :            : 
    2907         [ +  + ]:     828299 :   if (!signe(x))
    2908                 :            :   {
    2909                 :      24109 :     long e = expo(x);
    2910                 :      24109 :     *s = real_0_bit(e);
    2911         [ -  + ]:      24109 :     *c = e >= 0? real_0_bit(e): real_1(nbits2prec(-e));
    2912                 :     828299 :     return;
    2913                 :            :   }
    2914                 :            : 
    2915                 :     804190 :   av=avma; p1=mpcosm1(x,&mod8); tetpil=avma;
    2916   [ +  +  +  +  :     804190 :   switch(mod8)
             +  +  +  +  
                      - ]
    2917                 :            :   {
    2918                 :     218580 :     case 0: *c=addsr( 1,p1); *s=mpaut(p1); break;
    2919                 :      16932 :     case 1: *s=addsr( 1,p1); *c=mpaut(p1); togglesign(*c); break;
    2920                 :     204809 :     case 2: *c=subsr(-1,p1); *s=mpaut(p1); togglesign(*s); break;
    2921                 :      25139 :     case 3: *s=subsr(-1,p1); *c=mpaut(p1); break;
    2922                 :     170384 :     case 4: *c=addsr( 1,p1); *s=mpaut(p1); togglesign(*s); break;
    2923                 :      20566 :     case 5: *s=addsr( 1,p1); *c=mpaut(p1); break;
    2924                 :     104364 :     case 6: *c=subsr(-1,p1); *s=mpaut(p1); break;
    2925                 :      43416 :     case 7: *s=subsr(-1,p1); *c=mpaut(p1); togglesign(*c); break;
    2926                 :            :   }
    2927                 :     804190 :   gptr[0]=s; gptr[1]=c;
    2928                 :     804190 :   gerepilemanysp(av,tetpil,gptr,2);
    2929                 :            : }
    2930                 :            : 
    2931                 :            : /* SINE and COSINE - 1 */
    2932                 :            : void
    2933                 :         56 : mpsincosm1(GEN x, GEN *s, GEN *c)
    2934                 :            : {
    2935                 :            :   long mod8;
    2936                 :            :   pari_sp av, tetpil;
    2937                 :            :   GEN p1, *gptr[2];
    2938                 :            : 
    2939         [ -  + ]:         56 :   if (!signe(x))
    2940                 :            :   {
    2941                 :          0 :     long e = expo(x);
    2942                 :          0 :     *s = real_0_bit(e);
    2943                 :          0 :     *c = real_0_bit(2*e-1);
    2944                 :         56 :     return;
    2945                 :            :   }
    2946                 :         56 :   av=avma; p1=mpcosm1(x,&mod8); tetpil=avma;
    2947   [ -  -  -  -  :         56 :   switch(mod8)
             -  +  -  -  
                      - ]
    2948                 :            :   {
    2949                 :          0 :     case 0: *c=rcopy(p1); *s=mpaut(p1); break;
    2950                 :          0 :     case 1: *s=addsr(1,p1); *c=subrs(mpaut(p1),1); togglesign(*c); break;
    2951                 :          0 :     case 2: *c=subsr(-2,p1); *s=mpaut(p1); togglesign(*s); break;
    2952                 :          0 :     case 3: *s=subsr(-1,p1); *c=subrs(mpaut(p1),1); break;
    2953                 :          0 :     case 4: *c=rcopy(p1); *s=mpaut(p1); togglesign(*s); break;
    2954                 :         56 :     case 5: *s=addsr( 1,p1); *c=subrs(mpaut(p1),1); break;
    2955                 :          0 :     case 6: *c=subsr(-2,p1); *s=mpaut(p1); break;
    2956                 :          0 :     case 7: *s=subsr(-1,p1); *c=subsr(-1,mpaut(p1)); break;
    2957                 :            :   }
    2958                 :         56 :   gptr[0]=s; gptr[1]=c;
    2959                 :         56 :   gerepilemanysp(av,tetpil,gptr,2);
    2960                 :            : }
    2961                 :            : 
    2962                 :            : /* return exp(ix), x a t_REAL */
    2963                 :            : GEN
    2964                 :     102650 : expIr(GEN x)
    2965                 :            : {
    2966                 :     102650 :   pari_sp av = avma;
    2967                 :     102650 :   GEN v = cgetg(3,t_COMPLEX);
    2968                 :     102650 :   mpsincos(x, (GEN*)(v+2), (GEN*)(v+1));
    2969         [ +  + ]:     102650 :   if (!signe(gel(v,2))) return gerepilecopy(av, gel(v,1));
    2970                 :     102650 :   return v;
    2971                 :            : }
    2972                 :            : 
    2973                 :            : /* return exp(ix)-1, x a t_REAL */
    2974                 :            : static GEN
    2975                 :         56 : expm1_Ir(GEN x)
    2976                 :            : {
    2977                 :         56 :   pari_sp av = avma;
    2978                 :         56 :   GEN v = cgetg(3,t_COMPLEX);
    2979                 :         56 :   mpsincosm1(x, (GEN*)(v+2), (GEN*)(v+1));
    2980         [ -  + ]:         56 :   if (!signe(gel(v,2))) return gerepilecopy(av, gel(v,1));
    2981                 :         56 :   return v;
    2982                 :            : }
    2983                 :            : 
    2984                 :            : /* return exp(z)-1, z complex */
    2985                 :            : GEN
    2986                 :         63 : cxexpm1(GEN z, long prec)
    2987                 :            : {
    2988                 :         63 :   pari_sp av = avma;
    2989                 :         63 :   GEN X, Y, x = real_i(z), y = imag_i(z);
    2990         [ -  + ]:         63 :   if (typ(x) != t_REAL) x = gtofp(x, prec);
    2991         [ +  - ]:         63 :   if (typ(y) != t_REAL) y = gtofp(y, prec);
    2992         [ +  + ]:         63 :   if (gequal0(y)) return mpexpm1(x);
    2993         [ -  + ]:         56 :   if (gequal0(x)) return expm1_Ir(y);
    2994                 :         56 :   X = mpexpm1(x); /* t_REAL */
    2995                 :         56 :   Y = expm1_Ir(y);
    2996                 :            :   /* exp(x+iy) - 1 = (exp(x)-1)(exp(iy)-1) + exp(x)-1 + exp(iy)-1 */
    2997                 :         63 :   return gerepileupto(av, gadd(gadd(X,Y), gmul(X,Y)));
    2998                 :            : }
    2999                 :            : 
    3000                 :            : void
    3001                 :     716612 : gsincos(GEN x, GEN *s, GEN *c, long prec)
    3002                 :            : {
    3003                 :            :   long i, j, ex, ex2, lx, ly, mi;
    3004                 :            :   pari_sp av, tetpil;
    3005                 :            :   GEN y, r, u, v, u1, v1, p1, p2, p3, p4, ps, pc;
    3006                 :            :   GEN *gptr[4];
    3007                 :            : 
    3008   [ +  +  +  -  :     716612 :   switch(typ(x))
                      + ]
    3009                 :            :   {
    3010                 :            :     case t_INT: case t_FRAC:
    3011                 :        168 :       *s = cgetr(prec);
    3012                 :        168 :       *c = cgetr(prec); av = avma;
    3013                 :        168 :       mpsincos(tofp_safe(x, prec), &ps, &pc);
    3014                 :        168 :       affrr_fixlg(ps,*s);
    3015                 :        168 :       affrr_fixlg(pc,*c); avma = av; return;
    3016                 :            : 
    3017                 :            :     case t_REAL:
    3018                 :     716122 :       mpsincos(x,s,c); return;
    3019                 :            : 
    3020                 :            :     case t_COMPLEX:
    3021         [ +  + ]:        147 :       i = precision(x); if (!i) i = prec;
    3022                 :        147 :       ps = cgetc(i); *s = ps;
    3023                 :        147 :       pc = cgetc(i); *c = pc; av = avma;
    3024                 :        147 :       r = gexp(gel(x,2),prec);
    3025                 :        147 :       v1 = gmul2n(addrr(invr(r),r), -1); /* = cos(I*Im(x)) */
    3026                 :        147 :       u1 = subrr(r, v1); /* = I*sin(I*Im(x)) */
    3027                 :        147 :       gsincos(gel(x,1), &u,&v, prec);
    3028                 :        147 :       affrr_fixlg(mulrr(v1,u), gel(ps,1));
    3029                 :        147 :       affrr_fixlg(mulrr(u1,v), gel(ps,2));
    3030                 :        147 :       affrr_fixlg(mulrr(v1,v), gel(pc,1));
    3031                 :        147 :       affrr_fixlg(mulrr(u1,u), gel(pc,2)); togglesign(gel(pc,2));
    3032                 :        147 :       avma = av; return;
    3033                 :            : 
    3034                 :            :     case t_QUAD:
    3035                 :          0 :       av = avma; gsincos(quadtofp(x, prec), s, c, prec);
    3036                 :          0 :       gerepileall(av, 2, s, c); return;
    3037                 :            : 
    3038                 :            :     default:
    3039         [ -  + ]:        175 :       av = avma; if (!(y = toser_i(x))) break;
    3040         [ -  + ]:        175 :       if (gequal0(y)) { *s = gerepilecopy(av,y); *c = gaddsg(1,*s); return; }
    3041                 :            : 
    3042                 :        175 :       ex = valp(y); lx = lg(y); ex2 = 2*ex+2;
    3043         [ -  + ]:        175 :       if (ex < 0) pari_err_DOMAIN("gsincos","valuation", "<", gen_0, x);
    3044         [ +  + ]:        175 :       if (ex2 > lx)
    3045                 :            :       {
    3046         [ +  - ]:          7 :         *s = x == y? gcopy(y): gerepilecopy(av, y); av = avma;
    3047                 :          7 :         *c = gerepileupto(av, gsubsg(1, gdivgs(gsqr(y),2)));
    3048                 :            :         return;
    3049                 :            :       }
    3050         [ +  + ]:        168 :       if (!ex)
    3051                 :            :       {
    3052                 :         21 :         gsincos(serchop0(y),&u,&v,prec);
    3053                 :         21 :         gsincos(gel(y,2),&u1,&v1,prec);
    3054                 :         21 :         p1 = gmul(v1,v);
    3055                 :         21 :         p2 = gmul(u1,u);
    3056                 :         21 :         p3 = gmul(v1,u);
    3057                 :         21 :         p4 = gmul(u1,v); tetpil = avma;
    3058                 :         21 :         *c = gsub(p1,p2);
    3059                 :         21 :         *s = gadd(p3,p4);
    3060                 :         21 :         gptr[0]=s; gptr[1]=c;
    3061                 :         21 :         gerepilemanysp(av,tetpil,gptr,2);
    3062                 :            :         return;
    3063                 :            :       }
    3064                 :            : 
    3065                 :        147 :       ly = lx+2*ex;
    3066 [ +  + ][ +  - ]:       2009 :       mi = lx-1; while (mi>=3 && isrationalzero(gel(y,mi))) mi--;
    3067                 :        147 :       mi += ex-2;
    3068                 :        147 :       pc = cgetg(ly,t_SER); *c = pc;
    3069                 :        147 :       ps = cgetg(lx,t_SER); *s = ps;
    3070                 :        147 :       pc[1] = evalsigne(1) | _evalvalp(0) | evalvarn(varn(y));
    3071                 :        147 :       gel(pc,2) = gen_1; ps[1] = y[1];
    3072         [ +  + ]:        301 :       for (i=2; i<ex+2; i++) gel(ps,i) = gcopy(gel(y,i));
    3073         [ +  + ]:        308 :       for (i=3; i< ex2; i++) gel(pc,i) = gen_0;
    3074         [ +  + ]:       2156 :       for (i=ex2; i<ly; i++)
    3075                 :            :       {
    3076                 :       2009 :         long ii = i-ex;
    3077                 :       2009 :         av = avma; p1 = gen_0;
    3078         [ +  + ]:       4018 :         for (j=ex; j<=minss(ii-2,mi); j++)
    3079                 :       2009 :           p1 = gadd(p1, gmulgs(gmul(gel(y,j-ex+2),gel(ps,ii-j)),j));
    3080                 :       2009 :         gel(pc,i) = gerepileupto(av, gdivgs(p1,2-i));
    3081         [ +  + ]:       2009 :         if (ii < lx)
    3082                 :            :         {
    3083                 :       1855 :           av = avma; p1 = gen_0;
    3084         [ +  + ]:       3556 :           for (j=ex; j<=minss(i-ex2,mi); j++)
    3085                 :       1701 :             p1 = gadd(p1,gmulgs(gmul(gel(y,j-ex+2),gel(pc,i-j)),j));
    3086                 :       1855 :           p1 = gdivgs(p1,i-2);
    3087                 :       1855 :           gel(ps,ii) = gerepileupto(av, gadd(p1,gel(y,ii)));
    3088                 :            :         }
    3089                 :            :       }
    3090                 :            :       return;
    3091                 :            :   }
    3092                 :     716612 :   pari_err_TYPE("gsincos",x);
    3093                 :            : }
    3094                 :            : 
    3095                 :            : /********************************************************************/
    3096                 :            : /**                                                                **/
    3097                 :            : /**                     TANGENT and COTANGENT                      **/
    3098                 :            : /**                                                                **/
    3099                 :            : /********************************************************************/
    3100                 :            : static GEN
    3101                 :         21 : mptan(GEN x)
    3102                 :            : {
    3103                 :         21 :   pari_sp av = avma;
    3104                 :            :   GEN s, c;
    3105                 :            : 
    3106                 :         21 :   mpsincos(x,&s,&c);
    3107         [ -  + ]:         21 :   if (!signe(c))
    3108                 :          0 :     pari_err_DOMAIN("tan", "argument", "=", strtoGENstr("Pi/2 + kPi"),x);
    3109                 :         21 :   return gerepileuptoleaf(av, divrr(s,c));
    3110                 :            : }
    3111                 :            : 
    3112                 :            : GEN
    3113                 :         77 : gtan(GEN x, long prec)
    3114                 :            : {
    3115                 :            :   pari_sp av;
    3116                 :            :   GEN y, s, c;
    3117                 :            : 
    3118   [ +  +  +  +  :         77 :   switch(typ(x))
                      + ]
    3119                 :            :   {
    3120                 :         14 :     case t_REAL: return mptan(x);
    3121                 :            : 
    3122                 :            :     case t_COMPLEX: {
    3123         [ +  + ]:         14 :       if (isintzero(gel(x,1))) retmkcomplex(gen_0,gtanh(gel(x,2),prec));
    3124                 :          7 :       av = avma; y = mulcxmI(gtanh(mulcxI(x), prec)); /* tan x = -I th(I x) */
    3125                 :          7 :       gel(y,1) = gcopy(gel(y,1));
    3126                 :          7 :       return gerepileupto(av, y);
    3127                 :            :     }
    3128                 :            :     case t_INT: case t_FRAC:
    3129                 :          7 :       y = cgetr(prec); av = avma;
    3130                 :          7 :       affrr_fixlg(mptan(tofp_safe(x,prec)), y); avma = av; return y;
    3131                 :            : 
    3132                 :            :     case t_PADIC:
    3133                 :         14 :       av = avma;
    3134                 :         14 :       return gerepileupto(av, gdiv(gsin(x,prec), gcos(x,prec)));
    3135                 :            : 
    3136                 :            :     default:
    3137         [ +  + ]:         28 :       av = avma; if (!(y = toser_i(x))) break;
    3138         [ -  + ]:         21 :       if (gequal0(y)) return gerepilecopy(av, y);
    3139         [ +  + ]:         21 :       if (valp(y) < 0)
    3140                 :          7 :         pari_err_DOMAIN("tan","valuation", "<", gen_0, x);
    3141                 :         14 :       gsincos(y,&s,&c,prec);
    3142                 :         14 :       return gerepileupto(av, gdiv(s,c));
    3143                 :            :   }
    3144                 :         70 :   return trans_eval("tan",gtan,x,prec);
    3145                 :            : }
    3146                 :            : 
    3147                 :            : static GEN
    3148                 :         21 : mpcotan(GEN x)
    3149                 :            : {
    3150                 :         21 :   pari_sp av=avma, tetpil;
    3151                 :            :   GEN s,c;
    3152                 :            : 
    3153                 :         21 :   mpsincos(x,&s,&c); tetpil=avma;
    3154                 :         21 :   return gerepile(av,tetpil,divrr(c,s));
    3155                 :            : }
    3156                 :            : 
    3157                 :            : GEN
    3158                 :        105 : gcotan(GEN x, long prec)
    3159                 :            : {
    3160                 :            :   pari_sp av;
    3161                 :            :   GEN y, s, c;
    3162                 :            : 
    3163   [ +  +  +  +  :        105 :   switch(typ(x))
                      + ]
    3164                 :            :   {
    3165                 :            :     case t_REAL:
    3166                 :         14 :       return mpcotan(x);
    3167                 :            : 
    3168                 :            :     case t_COMPLEX:
    3169         [ +  + ]:         28 :       if (isintzero(gel(x,1))) {
    3170                 :         14 :         GEN z = cgetg(3, t_COMPLEX);
    3171                 :         14 :         gel(z,1) = gen_0;
    3172                 :         14 :         av = avma;
    3173                 :         14 :         gel(z,2) = gerepileupto(av, gneg(ginv(gtanh(gel(x,2),prec))));
    3174                 :         14 :         return z;
    3175                 :            :       }
    3176                 :         14 :       av = avma;
    3177                 :         14 :       gsincos(x,&s,&c,prec);
    3178                 :         14 :       return gerepileupto(av, gdiv(c,s));
    3179                 :            : 
    3180                 :            :     case t_INT: case t_FRAC:
    3181                 :          7 :       y = cgetr(prec); av = avma;
    3182                 :          7 :       affrr_fixlg(mpcotan(tofp_safe(x,prec)), y); avma = av; return y;
    3183                 :            : 
    3184                 :            :     case t_PADIC:
    3185                 :         14 :       av = avma;
    3186                 :         14 :       return gerepileupto(av, gdiv(gcos(x,prec), gsin(x,prec)));
    3187                 :            : 
    3188                 :            :     default:
    3189         [ +  + ]:         42 :       av = avma; if (!(y = toser_i(x))) break;
    3190         [ -  + ]:         35 :       if (gequal0(y)) pari_err_DOMAIN("cotan", "argument", "=", gen_0, y);
    3191         [ +  + ]:         35 :       if (valp(y) < 0) pari_err_DOMAIN("cotan","valuation", "<", gen_0, x);
    3192                 :         28 :       gsincos(y,&s,&c,prec);
    3193                 :         28 :       return gerepileupto(av, gdiv(c,s));
    3194                 :            :   }
    3195                 :         91 :   return trans_eval("cotan",gcotan,x,prec);
    3196                 :            : }

Generated by: LCOV version 1.9