Code coverage tests

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

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

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

LCOV - code coverage report
Current view: top level - basemath - trans2.c (source / functions) Hit Total Coverage
Test: PARI/GP v2.10.0 lcov report (development 22303-eb3e11d) Lines: 1009 1048 96.3 %
Date: 2018-04-21 06:16:28 Functions: 65 65 100.0 %
Legend: Lines: hit not hit

          Line data    Source code
       1             : /* Copyright (C) 2000  The PARI group.
       2             : 
       3             : This file is part of the PARI/GP package.
       4             : 
       5             : PARI/GP is free software; you can redistribute it and/or modify it under the
       6             : terms of the GNU General Public License as published by the Free Software
       7             : Foundation. It is distributed in the hope that it will be useful, but WITHOUT
       8             : ANY WARRANTY WHATSOEVER.
       9             : 
      10             : Check the License for details. You should have received a copy of it, along
      11             : with the package; see the file 'COPYING'. If not, write to the Free Software
      12             : Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. */
      13             : 
      14             : /********************************************************************/
      15             : /**                                                                **/
      16             : /**                   TRANSCENDENTAL FUNCTIONS                     **/
      17             : /**                          (part 2)                              **/
      18             : /**                                                                **/
      19             : /********************************************************************/
      20             : #include "pari.h"
      21             : #include "paripriv.h"
      22             : 
      23             : GEN
      24      188109 : trans_fix_arg(long *prec, GEN *s0, GEN *sig, GEN *tau, pari_sp *av, GEN *res)
      25             : {
      26             :   GEN s, p1;
      27             :   long l;
      28      188109 :   if (typ(*s0)==t_COMPLEX && gequal0(gel(*s0,2))) *s0 = gel(*s0,1);
      29      188109 :   s = *s0;
      30      188109 :   l = precision(s); if (!l) l = *prec;
      31      188109 :   if (l < LOWDEFAULTPREC) l = LOWDEFAULTPREC;
      32      188109 :   *res = cgetc(l); *av = avma;
      33      188109 :   if (typ(s) == t_COMPLEX)
      34             :   { /* s = sig + i t */
      35      175567 :     s = cxtofp(s, l+EXTRAPRECWORD);
      36      175567 :     *sig = gel(s,1);
      37      175567 :     *tau = gel(s,2);
      38             :   }
      39             :   else /* real number */
      40             :   {
      41       12542 :     *sig = s = gtofp(s, l+EXTRAPRECWORD);
      42       12542 :     *tau = gen_0;
      43       12542 :     p1 = trunc2nr(s, 0);
      44       12542 :     if (!signe(subri(s,p1))) *s0 = p1;
      45             :   }
      46      188109 :   *prec = l; return s;
      47             : }
      48             : 
      49             : /********************************************************************/
      50             : /**                                                                **/
      51             : /**                          ARCTANGENT                            **/
      52             : /**                                                                **/
      53             : /********************************************************************/
      54             : /* atan(b/a), real a and b, suitable for gerepileupto */
      55             : static GEN
      56         196 : atan2_agm(GEN a, GEN b, long prec)
      57         196 : { return gel(logagmcx(mkcomplex(a, b), prec), 2); }
      58             : static GEN
      59      504870 : mpatan(GEN x)
      60             : {
      61      504870 :   long l, l1, l2, n, m, i, lp, e, s, sx = signe(x);
      62             :   pari_sp av0, av;
      63             :   double alpha, beta, delta;
      64             :   GEN y, p1, p2, p3, p4, p5, unr;
      65             :   int inv;
      66             : 
      67      504870 :   if (!sx) return real_0_bit(expo(x));
      68      504807 :   l = lp = realprec(x);
      69      504807 :   if (absrnz_equal1(x)) { /* |x| = 1 */
      70        3677 :     y = Pi2n(-2, l+EXTRAPRECWORD); if (sx < 0) setsigne(y,-1);
      71        3677 :     return y;
      72             :   }
      73      501130 :   if (l > AGM_ATAN_LIMIT)
      74         175 :   { av = avma; return gerepileuptoleaf(av, atan2_agm(gen_1, x, l)); }
      75             : 
      76      500955 :   e = expo(x); inv = (e >= 0); /* = (|x| > 1 ) */
      77      500955 :   if (e > 0) lp += nbits2extraprec(e);
      78             : 
      79      500955 :   y = cgetr(lp); av0 = avma;
      80      500955 :   p1 = rtor(x, l+EXTRAPRECWORD); setabssign(p1); /* p1 = |x| */
      81      500955 :   if (inv) p1 = invr(p1);
      82      500955 :   e = expo(p1);
      83      500955 :   if (e < -100)
      84        2760 :     alpha = 1.65149612947 - e; /* log_2(Pi) - e */
      85             :   else
      86      498195 :     alpha = log2(M_PI / atan(rtodbl(p1)));
      87      500955 :   beta = (double)(prec2nbits(l)>>1);
      88      500955 :   delta = 1 + beta - alpha/2;
      89      500955 :   if (delta <= 0) { n = 1; m = 0; }
      90             :   else
      91             :   {
      92      498240 :     double fi = alpha-2;
      93      498240 :     if (delta >= fi*fi)
      94             :     {
      95      489484 :       double t = 1 + sqrt(delta);
      96      489484 :       n = (long)t;
      97      489484 :       m = (long)(t - fi);
      98             :     }
      99             :     else
     100             :     {
     101        8756 :       n = (long)(1+beta/fi);
     102        8756 :       m = 0;
     103             :     }
     104             :   }
     105      500955 :   l2 = l + nbits2extraprec(m);
     106      500955 :   p2 = rtor(p1, l2); av = avma;
     107     4929154 :   for (i=1; i<=m; i++)
     108             :   {
     109     4428199 :     p5 = addsr(1, sqrr(p2)); setprec(p5,l2);
     110     4428199 :     p5 = addsr(1, sqrtr_abs(p5)); setprec(p5,l2);
     111     4428199 :     affrr(divrr(p2,p5), p2); avma = av;
     112             :   }
     113      500955 :   p3 = sqrr(p2); l1 = minss(LOWDEFAULTPREC+EXTRAPRECWORD, l2); /* l1 increases to l2 */;
     114      500955 :   unr = real_1(l2); setprec(unr,l1);
     115      500955 :   p4 = cgetr(l2); setprec(p4,l1);
     116      500955 :   affrr(divru(unr,2*n+1), p4);
     117      500955 :   s = 0; e = expo(p3); av = avma;
     118     5406838 :   for (i = n; i > 1; i--) /* n >= 1. i = 1 done outside for efficiency */
     119             :   {
     120     4905883 :     setprec(p3,l1); p5 = mulrr(p4,p3);
     121     4905883 :     l1 += dvmdsBIL(s - e, &s); if (l1 > l2) l1 = l2;
     122     4905883 :     setprec(unr,l1); p5 = subrr(divru(unr,2*i-1), p5);
     123     4905883 :     setprec(p4,l1); affrr(p5,p4); avma = av;
     124             :   }
     125      500955 :   setprec(p3, l2); p5 = mulrr(p4,p3); /* i = 1 */
     126      500955 :   setprec(unr,l2); p4 = subrr(unr, p5);
     127             : 
     128      500955 :   p4 = mulrr(p2,p4); shiftr_inplace(p4, m);
     129      500955 :   if (inv) p4 = subrr(Pi2n(-1, lp), p4);
     130      500955 :   if (sx < 0) togglesign(p4);
     131      500955 :   affrr_fixlg(p4,y); avma = av0; return y;
     132             : }
     133             : 
     134             : GEN
     135       19318 : gatan(GEN x, long prec)
     136             : {
     137             :   pari_sp av;
     138             :   GEN a, y;
     139             : 
     140       19318 :   switch(typ(x))
     141             :   {
     142       11086 :     case t_REAL: return mpatan(x);
     143             :     case t_COMPLEX: /* atan(x) = -i atanh(ix) */
     144        7672 :       if (ismpzero(gel(x,2))) return gatan(gel(x,1), prec);
     145        7441 :       av = avma; return gerepilecopy(av, mulcxmI(gatanh(mulcxI(x),prec)));
     146             :     default:
     147         560 :       av = avma; if (!(y = toser_i(x))) break;
     148          28 :       if (valp(y) < 0) pari_err_DOMAIN("atan","valuation", "<", gen_0, x);
     149          21 :       if (lg(y)==2) return gerepilecopy(av, y);
     150             :       /* lg(y) > 2 */
     151          14 :       a = integser(gdiv(derivser(y), gaddsg(1,gsqr(y))));
     152          14 :       if (!valp(y)) a = gadd(a, gatan(gel(y,2),prec));
     153          14 :       return gerepileupto(av, a);
     154             :   }
     155         532 :   return trans_eval("atan",gatan,x,prec);
     156             : }
     157             : /********************************************************************/
     158             : /**                                                                **/
     159             : /**                             ARCSINE                            **/
     160             : /**                                                                **/
     161             : /********************************************************************/
     162             : /* |x| < 1, x != 0 */
     163             : static GEN
     164          98 : mpasin(GEN x) {
     165          98 :   pari_sp av = avma;
     166          98 :   GEN z, a = sqrtr(subsr(1, sqrr(x)));
     167          98 :   if (realprec(x) > AGM_ATAN_LIMIT)
     168           7 :     z = atan2_agm(a, x, realprec(x));
     169             :   else
     170          91 :     z = mpatan(divrr(x, a));
     171          98 :   return gerepileuptoleaf(av, z);
     172             : }
     173             : 
     174             : static GEN mpacosh(GEN x);
     175             : GEN
     176        8183 : gasin(GEN x, long prec)
     177             : {
     178             :   long sx;
     179             :   pari_sp av;
     180             :   GEN a, y, p1;
     181             : 
     182        8183 :   switch(typ(x))
     183             :   {
     184         483 :     case t_REAL: sx = signe(x);
     185         483 :       if (!sx) return real_0_bit(expo(x));
     186         476 :       if (absrnz_equal1(x)) { /* |x| = 1 */
     187          28 :         if (sx > 0) return Pi2n(-1, realprec(x)); /* 1 */
     188          14 :         y = Pi2n(-1, realprec(x)); setsigne(y, -1); return y; /* -1 */
     189             :       }
     190         448 :       if (expo(x) < 0) return mpasin(x);
     191         350 :       y = cgetg(3,t_COMPLEX);
     192         350 :       gel(y,1) = Pi2n(-1, realprec(x));
     193         350 :       gel(y,2) = mpacosh(x);
     194         350 :       if (sx < 0) togglesign(gel(y,1)); else togglesign(gel(y,2));
     195         350 :       return y;
     196             : 
     197             :     case t_COMPLEX: /* asin(z) = -i asinh(iz) */
     198        7637 :       if (ismpzero(gel(x,2))) return gasin(gel(x,1), prec);
     199        7406 :       av = avma;
     200        7406 :       return gerepilecopy(av, mulcxmI(gasinh(mulcxI(x), prec)));
     201             :     default:
     202          63 :       av = avma; if (!(y = toser_i(x))) break;
     203          42 :       if (gequal0(y)) return gerepilecopy(av, y);
     204             :       /* lg(y) > 2*/
     205          35 :       if (valp(y) < 0) pari_err_DOMAIN("asin","valuation", "<", gen_0, x);
     206          28 :       p1 = gsubsg(1,gsqr(y));
     207          28 :       if (gequal0(p1))
     208             :       {
     209          21 :         GEN t = Pi2n(-1,prec);
     210          21 :         if (gsigne(gel(y,2)) < 0) setsigne(t, -1);
     211          21 :         return gerepileupto(av, scalarser(t, varn(y), valp(p1)>>1));
     212             :       }
     213           7 :       p1 = gdiv(derivser(y), gsqrt(p1,prec));
     214           7 :       a = integser(p1);
     215           7 :       if (!valp(y)) a = gadd(a, gasin(gel(y,2),prec));
     216           7 :       return gerepileupto(av, a);
     217             :   }
     218          21 :   return trans_eval("asin",gasin,x,prec);
     219             : }
     220             : /********************************************************************/
     221             : /**                                                                **/
     222             : /**                             ARCCOSINE                          **/
     223             : /**                                                                **/
     224             : /********************************************************************/
     225             : static GEN
     226          14 : acos0(long e) { return Pi2n(-1, nbits2prec(e<0? -e: 1)); }
     227             : 
     228             : /* |x| < 1, x != 0 */
     229             : static GEN
     230         105 : mpacos(GEN x)
     231             : {
     232         105 :   pari_sp av = avma;
     233         105 :   GEN z, a = sqrtr(subsr(1, sqrr(x)));
     234         105 :   if (realprec(x) > AGM_ATAN_LIMIT)
     235          14 :     z = atan2_agm(x, a, realprec(x));
     236             :   else
     237             :   {
     238          91 :     z = mpatan(divrr(a, x));
     239          91 :     if (signe(x) < 0) z = addrr(mppi(realprec(z)), z);
     240             :   }
     241         105 :   return gerepileuptoleaf(av, z);
     242             : }
     243             : 
     244             : GEN
     245        7938 : gacos(GEN x, long prec)
     246             : {
     247             :   long sx;
     248             :   pari_sp av;
     249             :   GEN a, y, p1;
     250             : 
     251        7938 :   switch(typ(x))
     252             :   {
     253         252 :     case t_REAL: sx = signe(x);
     254         252 :       if (!sx) return acos0(expo(x));
     255         245 :       if (absrnz_equal1(x)) /* |x| = 1 */
     256          14 :         return sx > 0? real_0_bit( -(bit_prec(x)>>1) ) : mppi(realprec(x));
     257         231 :       if (expo(x) < 0) return mpacos(x);
     258             : 
     259         175 :       y = cgetg(3,t_COMPLEX); p1 = mpacosh(x);
     260         175 :       if (sx < 0) { gel(y,1) = mppi(realprec(x)); togglesign(p1); }
     261          91 :       else gel(y,1) = gen_0;
     262         175 :       gel(y,2) = p1; return y;
     263             : 
     264             :     case t_COMPLEX:
     265        7637 :       if (ismpzero(gel(x,2))) return gacos(gel(x,1), prec);
     266        7406 :       av = avma;
     267        7406 :       p1 = gadd(x, mulcxI(gsqrt(gsubsg(1,gsqr(x)), prec)));
     268        7406 :       y = glog(p1,prec); /* log(x + I*sqrt(1-x^2)) */
     269        7406 :       return gerepilecopy(av, mulcxmI(y));
     270             :     default:
     271          49 :       av = avma; if (!(y = toser_i(x))) break;
     272          35 :       if (valp(y) < 0) pari_err_DOMAIN("acos","valuation", "<", gen_0, x);
     273          28 :       if (lg(y) > 2)
     274             :       {
     275          21 :         p1 = gsubsg(1,gsqr(y));
     276          21 :         if (gequal0(p1)) { avma = av; return zeroser(varn(y), valp(p1)>>1); }
     277           7 :         p1 = integser(gdiv(gneg(derivser(y)), gsqrt(p1,prec)));
     278             :         /*y(t) = 1+O(t)*/
     279           7 :         if (gequal1(gel(y,2)) && !valp(y)) return gerepileupto(av, p1);
     280             :       }
     281           7 :       else p1 = y;
     282          14 :       a = (lg(y)==2 || valp(y))? Pi2n(-1, prec): gacos(gel(y,2),prec);
     283          14 :       return gerepileupto(av, gadd(a,p1));
     284             :   }
     285          14 :   return trans_eval("acos",gacos,x,prec);
     286             : }
     287             : /********************************************************************/
     288             : /**                                                                **/
     289             : /**                            ARGUMENT                            **/
     290             : /**                                                                **/
     291             : /********************************************************************/
     292             : 
     293             : /* we know that x and y are not both 0 */
     294             : static GEN
     295      493728 : mparg(GEN x, GEN y)
     296             : {
     297      493728 :   long prec, sx = signe(x), sy = signe(y);
     298             :   GEN z;
     299             : 
     300      493728 :   if (!sy)
     301             :   {
     302         112 :     if (sx > 0) return real_0_bit(expo(y) - expo(x));
     303          49 :     return mppi(realprec(x));
     304             :   }
     305      493616 :   prec = realprec(y); if (prec < realprec(x)) prec = realprec(x);
     306      493616 :   if (!sx)
     307             :   {
     308          14 :     z = Pi2n(-1, prec); if (sy < 0) setsigne(z,-1);
     309          14 :     return z;
     310             :   }
     311             : 
     312      493602 :   if (expo(x)-expo(y) > -2)
     313             :   {
     314      431915 :     z = mpatan(divrr(y,x)); if (sx > 0) return z;
     315      105979 :     return addrr_sign(z, signe(z), mppi(prec), sy);
     316             :   }
     317       61687 :   z = mpatan(divrr(x,y));
     318       61687 :   return addrr_sign(z, -signe(z), Pi2n(-1, prec), sy);
     319             : }
     320             : 
     321             : static GEN
     322      987456 : rfix(GEN x,long prec)
     323             : {
     324      987456 :   switch(typ(x))
     325             :   {
     326        2803 :     case t_INT: return itor(x, prec);
     327        5285 :     case t_FRAC: return fractor(x, prec);
     328      979368 :     case t_REAL: break;
     329           0 :     default: pari_err_TYPE("rfix (conversion to t_REAL)",x);
     330             :   }
     331      979368 :   return x;
     332             : }
     333             : 
     334             : static GEN
     335      493728 : cxarg(GEN x, GEN y, long prec)
     336             : {
     337      493728 :   pari_sp av = avma;
     338      493728 :   x = rfix(x,prec);
     339      493728 :   y = rfix(y,prec); return gerepileuptoleaf(av, mparg(x,y));
     340             : }
     341             : 
     342             : GEN
     343      494869 : garg(GEN x, long prec)
     344             : {
     345             :   long l;
     346      494869 :   if (gequal0(x)) pari_err_DOMAIN("arg", "argument", "=", gen_0, x);
     347      494869 :   switch(typ(x))
     348             :   {
     349        1141 :     case t_REAL: prec = realprec(x); /* fall through */
     350        1141 :     case t_INT: case t_FRAC: return (gsigne(x)>0)? real_0(prec): mppi(prec);
     351             :     case t_COMPLEX:
     352      493728 :       l = precision(x); if (l) prec = l;
     353      493728 :       return cxarg(gel(x,1),gel(x,2),prec);
     354             :   }
     355           0 :   return trans_eval("arg",garg,x,prec);
     356             : }
     357             : 
     358             : /********************************************************************/
     359             : /**                                                                **/
     360             : /**                      HYPERBOLIC COSINE                         **/
     361             : /**                                                                **/
     362             : /********************************************************************/
     363             : 
     364             : static GEN
     365        3584 : mpcosh(GEN x)
     366             : {
     367             :   pari_sp av;
     368             :   GEN z;
     369             : 
     370        3584 :   if (!signe(x)) { /* 1 + x */
     371           7 :     long e = expo(x);
     372           7 :     return e >= 0? real_0_bit(e): real_1_bit(-e);
     373             :   }
     374        3577 :   av = avma;
     375        3577 :   z = mpexp(x); z = addrr(z, invr(z)); shiftr_inplace(z, -1);
     376        3577 :   return gerepileuptoleaf(av, z);
     377             : }
     378             : 
     379             : GEN
     380        3668 : gcosh(GEN x, long prec)
     381             : {
     382             :   pari_sp av;
     383             :   GEN y, p1;
     384             : 
     385        3668 :   switch(typ(x))
     386             :   {
     387        3584 :     case t_REAL: return mpcosh(x);
     388             :     case t_COMPLEX:
     389          21 :       if (isintzero(gel(x,1))) return gcos(gel(x,2),prec);
     390             :       /* fall through */
     391             :     case t_PADIC:
     392          21 :       av = avma; p1 = gexp(x,prec); p1 = gadd(p1, ginv(p1));
     393          21 :       return gerepileupto(av, gmul2n(p1,-1));
     394             :     default:
     395          49 :       av = avma; if (!(y = toser_i(x))) break;
     396          28 :       if (gequal0(y) && valp(y) == 0) return gerepilecopy(av, y);
     397          28 :       p1 = gexp(y,prec); p1 = gadd(p1, ginv(p1));
     398          28 :       return gerepileupto(av, gmul2n(p1,-1));
     399             :   }
     400          21 :   return trans_eval("cosh",gcosh,x,prec);
     401             : }
     402             : /********************************************************************/
     403             : /**                                                                **/
     404             : /**                       HYPERBOLIC SINE                          **/
     405             : /**                                                                **/
     406             : /********************************************************************/
     407             : 
     408             : static GEN
     409          70 : mpsinh(GEN x)
     410             : {
     411             :   pari_sp av;
     412          70 :   long ex = expo(x), lx;
     413             :   GEN z, res;
     414             : 
     415          70 :   if (!signe(x)) return real_0_bit(ex);
     416          70 :   lx = realprec(x); res = cgetr(lx); av = avma;
     417          70 :   if (ex < 1 - BITS_IN_LONG)
     418             :   { /* y = e^x-1; e^x - e^(-x) = y(1 + 1/(y+1)) */
     419           7 :     GEN y = mpexpm1(x);
     420           7 :     z = addrs(y,1); if (lg(z) > lx+1) z = rtor(z,lx+1); /* e^x */
     421           7 :     z = mulrr(y, addsr(1,invr(z)));
     422             :   }
     423             :   else
     424             :   {
     425          63 :     z = mpexp(x);
     426          63 :     z = subrr(z, invr(z));
     427             :   }
     428          70 :   shiftr_inplace(z, -1);
     429          70 :   affrr(z, res); avma = av; return res;
     430             : }
     431             : 
     432             : GEN
     433        1792 : gsinh(GEN x, long prec)
     434             : {
     435             :   pari_sp av;
     436             :   GEN y, p1;
     437             : 
     438        1792 :   switch(typ(x))
     439             :   {
     440          70 :     case t_REAL: return mpsinh(x);
     441             :     case t_COMPLEX:
     442          21 :       if (isintzero(gel(x,1))) retmkcomplex(gen_0, gsin(gel(x,2),prec));
     443             :       /* fall through */
     444             :     case t_PADIC:
     445          14 :       av = avma; p1 = gexp(x,prec); p1 = gsub(p1, ginv(p1));
     446          14 :       return gerepileupto(av, gmul2n(p1,-1));
     447             :     default:
     448        1694 :       av = avma; if (!(y = toser_i(x))) break;
     449        1666 :       if (gequal0(y) && valp(y) == 0) return gerepilecopy(av, y);
     450        1666 :       p1 = gexp(y, prec); p1 = gsub(p1, ginv(p1));
     451        1666 :       return gerepileupto(av, gmul2n(p1,-1));
     452             :   }
     453          28 :   return trans_eval("sinh",gsinh,x,prec);
     454             : }
     455             : /********************************************************************/
     456             : /**                                                                **/
     457             : /**                      HYPERBOLIC TANGENT                        **/
     458             : /**                                                                **/
     459             : /********************************************************************/
     460             : 
     461             : static GEN
     462       77056 : mptanh(GEN x)
     463             : {
     464       77056 :   long lx, s = signe(x);
     465             :   GEN y;
     466             : 
     467       77056 :   if (!s) return real_0_bit(expo(x));
     468       77056 :   lx = realprec(x);
     469       77056 :   if (abscmprr(x, stor(prec2nbits(lx), LOWDEFAULTPREC)) >= 0) {
     470       24840 :     y = real_1(lx);
     471             :   } else {
     472       52216 :     pari_sp av = avma;
     473       52216 :     long ex = expo(x);
     474             :     GEN t;
     475       52216 :     if (ex < 1 - BITS_IN_LONG) x = rtor(x, lx + nbits2extraprec(-ex)-1);
     476       52216 :     t = exp1r_abs(gmul2n(x,1)); /* exp(|2x|) - 1 */
     477       52216 :     y = gerepileuptoleaf(av, divrr(t, addsr(2,t)));
     478             :   }
     479       77056 :   if (s < 0) togglesign(y); /* tanh is odd */
     480       77056 :   return y;
     481             : }
     482             : 
     483             : GEN
     484       77161 : gtanh(GEN x, long prec)
     485             : {
     486             :   pari_sp av;
     487             :   GEN y, t;
     488             : 
     489       77161 :   switch(typ(x))
     490             :   {
     491       77056 :     case t_REAL: return mptanh(x);
     492             :     case t_COMPLEX:
     493          35 :       if (isintzero(gel(x,1))) retmkcomplex(gen_0, gtan(gel(x,2),prec));
     494             :       /* fall through */
     495             :     case t_PADIC:
     496          28 :       av = avma;
     497          28 :       t = gexp(gmul2n(x,1),prec);
     498          28 :       t = gdivsg(-2, gaddgs(t,1));
     499          28 :       return gerepileupto(av, gaddsg(1,t));
     500             :     default:
     501          63 :       av = avma; if (!(y = toser_i(x))) break;
     502          28 :       if (gequal0(y)) return gerepilecopy(av, y);
     503          14 :       t = gexp(gmul2n(y, 1),prec);
     504          14 :       t = gdivsg(-2, gaddgs(t,1));
     505          14 :       return gerepileupto(av, gaddsg(1,t));
     506             :   }
     507          35 :   return trans_eval("tanh",gtanh,x,prec);
     508             : }
     509             : 
     510             : static GEN
     511           7 : mpcotanh(GEN x)
     512             : {
     513           7 :   long lx, s = signe(x);
     514             :   GEN y;
     515             : 
     516           7 :   if (!s) pari_err_DOMAIN("cotan", "argument", "=", gen_0, x);
     517             : 
     518           7 :   lx = realprec(x);
     519           7 :   if (abscmprr(x, stor(prec2nbits(lx), LOWDEFAULTPREC)) >= 0) {
     520           0 :     y = real_1(lx);
     521             :   } else {
     522           7 :     pari_sp av = avma;
     523           7 :     long ex = expo(x);
     524             :     GEN t;
     525           7 :     if (ex < 1 - BITS_IN_LONG) x = rtor(x, lx + nbits2extraprec(-ex)-1);
     526           7 :     t = exp1r_abs(gmul2n(x,1)); /* exp(|2x|) - 1 */
     527           7 :     y = gerepileuptoleaf(av, divrr(addsr(2,t), t));
     528             :   }
     529           7 :   if (s < 0) togglesign(y); /* cotanh is odd */
     530           7 :   return y;
     531             : }
     532             : 
     533             : GEN
     534          63 : gcotanh(GEN x, long prec)
     535             : {
     536             :   pari_sp av;
     537             :   GEN y, t;
     538             : 
     539          63 :   switch(typ(x))
     540             :   {
     541           7 :     case t_REAL: return mpcotanh(x);
     542             :     case t_COMPLEX:
     543          14 :       if (isintzero(gel(x,1))) retmkcomplex(gen_0, gcotan(gel(x,2),prec));
     544             :       /* fall through */
     545             :     case t_PADIC:
     546          14 :       av = avma;
     547          14 :       t = gexpm1(gmul2n(x,1),prec);
     548          14 :       return gerepileupto(av, gaddsg(1, gdivsg(2,t)));
     549             :     default:
     550          35 :       av = avma; if (!(y = toser_i(x))) break;
     551          28 :       if (gequal0(y)) return gerepilecopy(av, y);
     552          14 :       t = gexpm1(gmul2n(y,1),prec);
     553          14 :       return gerepileupto(av, gaddsg(1, gdivsg(2,t)));
     554             :   }
     555           7 :   return trans_eval("cotanh",gcotanh,x,prec);
     556             : }
     557             : 
     558             : /********************************************************************/
     559             : /**                                                                **/
     560             : /**                     AREA HYPERBOLIC SINE                       **/
     561             : /**                                                                **/
     562             : /********************************************************************/
     563             : 
     564             : /* x != 0 */
     565             : static GEN
     566         483 : mpasinh(GEN x)
     567             : {
     568             :   GEN z, res;
     569             :   pari_sp av;
     570         483 :   long lx = realprec(x), ex = expo(x);
     571             : 
     572         483 :   res = cgetr(lx); av = avma;
     573         483 :   if (ex < 1 - BITS_IN_LONG) x = rtor(x, lx + nbits2extraprec(-ex)-1);
     574         483 :   z = logr_abs( addrr_sign(x,1, sqrtr_abs( addrs(sqrr(x), 1) ), 1) );
     575         483 :   if (signe(x) < 0) togglesign(z);
     576         483 :   affrr(z, res); avma = av; return res;
     577             : }
     578             : 
     579             : GEN
     580       15715 : gasinh(GEN x, long prec)
     581             : {
     582             :   pari_sp av;
     583             :   GEN a, y, p1;
     584             : 
     585       15715 :   switch(typ(x))
     586             :   {
     587             :     case t_REAL:
     588         490 :       if (!signe(x)) return rcopy(x);
     589         483 :       return mpasinh(x);
     590             : 
     591             :     case t_COMPLEX:
     592       15036 :       if (ismpzero(gel(x,2))) return gasinh(gel(x,1), prec);
     593       14581 :       av = avma;
     594       14581 :       if (ismpzero(gel(x,1))) /* avoid cancellation */
     595         231 :         return gerepilecopy(av, mulcxI(gasin(gel(x,2), prec)));
     596       14350 :       p1 = gadd(x, gsqrt(gaddsg(1,gsqr(x)), prec));
     597       14350 :       y = glog(p1,prec); /* log (x + sqrt(1+x^2)) */
     598       14350 :       return gerepileupto(av, y);
     599             :     default:
     600         189 :       av = avma; if (!(y = toser_i(x))) break;
     601         161 :       if (gequal0(y)) return gerepilecopy(av, y);
     602         154 :       if (valp(y) < 0) pari_err_DOMAIN("asinh","valuation", "<", gen_0, x);
     603         147 :       p1 = gaddsg(1,gsqr(y));
     604         147 :       if (gequal0(p1))
     605             :       {
     606          14 :         GEN t = PiI2n(-1,prec);
     607          14 :         if ( gsigne(imag_i(gel(y,2))) < 0 ) setsigne(gel(t,2), -1);
     608          14 :         return gerepileupto(av, scalarser(t, varn(y), valp(p1)>>1));
     609             :       }
     610         133 :       p1 = gdiv(derivser(y), gsqrt(p1,prec));
     611         133 :       a = integser(p1);
     612         133 :       if (!valp(y)) a = gadd(a, gasinh(gel(y,2),prec));
     613         133 :       return gerepileupto(av, a);
     614             :   }
     615          28 :   return trans_eval("asinh",gasinh,x,prec);
     616             : }
     617             : /********************************************************************/
     618             : /**                                                                **/
     619             : /**                     AREA HYPERBOLIC COSINE                     **/
     620             : /**                                                                **/
     621             : /********************************************************************/
     622             : 
     623             : /* |x| >= 1, return ach(|x|) */
     624             : static GEN
     625         728 : mpacosh(GEN x)
     626             : {
     627         728 :   pari_sp av = avma;
     628             :   GEN z;
     629         728 :   if (absrnz_equal1(x)) return real_0_bit(- bit_prec(x) >> 1);
     630         721 :   z = logr_abs( addrr_sign(x, 1, sqrtr( subrs(sqrr(x), 1) ), 1) );
     631         721 :   return gerepileuptoleaf(av, z);
     632             : }
     633             : 
     634             : GEN
     635        7980 : gacosh(GEN x, long prec)
     636             : {
     637             :   pari_sp av;
     638             :   GEN y, p1;
     639             : 
     640        7980 :   switch(typ(x))
     641             :   {
     642             :     case t_REAL: {
     643         266 :       long s = signe(x), e = expo(x);
     644             :       GEN a, b;
     645         266 :       if (s > 0 && e >= 0) return mpacosh(x);
     646             :       /* x < 1 */
     647         147 :       y = cgetg(3,t_COMPLEX); a = gen_0;
     648         147 :       if (s == 0) b = acos0(e);
     649         140 :       else if (e < 0) b = mpacos(x); /* -1 < x < 1 */
     650             :       else {
     651          91 :         if (!absrnz_equal1(x)) a = mpacosh(x);
     652          91 :         b = mppi(realprec(x));
     653             :       }
     654         147 :       gel(y,1) = a;
     655         147 :       gel(y,2) = b; return y;
     656             :     }
     657             :     case t_COMPLEX:
     658        7637 :       if (ismpzero(gel(x,2))) return gacosh(gel(x,1), prec);
     659        7406 :       av = avma;
     660        7406 :       p1 = gadd(x, gsqrt(gaddsg(-1,gsqr(x)), prec));
     661        7406 :       y = glog(p1,prec); /* log(x + sqrt(x^2-1)) */
     662        7406 :       if (signe(real_i(y)) < 0) y = gneg(y);
     663        7406 :       return gerepileupto(av, y);
     664             :     default: {
     665             :       GEN a;
     666             :       long v;
     667          77 :       av = avma; if (!(y = toser_i(x))) break;
     668          49 :       v = valp(y);
     669          49 :       if (v < 0) pari_err_DOMAIN("acosh","valuation", "<", gen_0, x);
     670          42 :       if (gequal0(y))
     671             :       {
     672           7 :         if (!v) return gerepilecopy(av, y);
     673           7 :         return gerepileupto(av, gadd(y, PiI2n(-1, prec)));
     674             :       }
     675          35 :       p1 = gsubgs(gsqr(y),1);
     676          35 :       if (gequal0(p1)) { avma = av; return zeroser(varn(y), valp(p1)>>1); }
     677          21 :       p1 = gdiv(derivser(y), gsqrt(p1,prec));
     678          21 :       a = integser(p1);
     679          21 :       if (v)
     680           7 :         p1 = PiI2n(-1, prec); /* I Pi/2 */
     681             :       else
     682             :       {
     683          14 :         p1 = gel(y,2); if (gequal1(p1)) return gerepileupto(av,a);
     684           7 :         p1 = gacosh(p1, prec);
     685             :       }
     686          14 :       return gerepileupto(av, gadd(p1,a));
     687             :     }
     688             :   }
     689          28 :   return trans_eval("acosh",gacosh,x,prec);
     690             : }
     691             : /********************************************************************/
     692             : /**                                                                **/
     693             : /**                     AREA HYPERBOLIC TANGENT                    **/
     694             : /**                                                                **/
     695             : /********************************************************************/
     696             : 
     697             : /* |x| < 1, x != 0 */
     698             : static GEN
     699          98 : mpatanh(GEN x)
     700             : {
     701          98 :   pari_sp av = avma;
     702          98 :   long ex = expo(x);
     703             :   GEN z;
     704          98 :   if (ex < 1 - BITS_IN_LONG) x = rtor(x, realprec(x) + nbits2extraprec(-ex)-1);
     705          98 :   z = invr( subsr(1,x) ); shiftr_inplace(z, 1); /* 2/(1-x)*/
     706          98 :   z = logr_abs( addrs(z,-1) );
     707          98 :   shiftr_inplace(z, -1); return gerepileuptoleaf(av, z);
     708             : }
     709             : 
     710             : GEN
     711       15603 : gatanh(GEN x, long prec)
     712             : {
     713             :   long sx;
     714             :   pari_sp av;
     715             :   GEN a, y, z;
     716             : 
     717       15603 :   switch(typ(x))
     718             :   {
     719             :     case t_REAL:
     720         483 :       sx = signe(x);
     721         483 :       if (!sx) return real_0_bit(expo(x));
     722         476 :       if (expo(x) < 0) return mpatanh(x);
     723             : 
     724         378 :       y = cgetg(3,t_COMPLEX);
     725         378 :       av = avma;
     726         378 :       z = subrs(x,1);
     727         378 :       if (!signe(z)) pari_err_DOMAIN("atanh", "argument", "=", gen_1, x);
     728         364 :       z = invr(z); shiftr_inplace(z, 1); /* 2/(x-1)*/
     729         364 :       z = addrs(z,1);
     730         364 :       if (!signe(z)) pari_err_DOMAIN("atanh", "argument", "=", gen_m1, x);
     731         350 :       z = logr_abs(z);
     732         350 :       shiftr_inplace(z, -1); /* (1/2)log((1+x)/(x-1)) */
     733         350 :       gel(y,1) = gerepileuptoleaf(av, z);
     734         350 :       gel(y,2) = Pi2n(-1, realprec(x));
     735         350 :       if (sx > 0) togglesign(gel(y,2));
     736         350 :       return y;
     737             : 
     738             :     case t_COMPLEX: /* 2/(1-z) - 1 = (1+z) / (1-z) */
     739       15071 :       if (ismpzero(gel(x,2))) return gatanh(gel(x,1), prec);
     740       14616 :       av = avma; z = glog( gaddgs(gdivsg(2,gsubsg(1,x)),-1), prec );
     741       14616 :       return gerepileupto(av, gmul2n(z,-1));
     742             : 
     743             :     default:
     744          49 :       av = avma; if (!(y = toser_i(x))) break;
     745          28 :       if (valp(y) < 0) pari_err_DOMAIN("atanh","valuation", "<", gen_0, x);
     746          21 :       z = gdiv(derivser(y), gsubsg(1,gsqr(y)));
     747          14 :       a = integser(z);
     748          14 :       if (!valp(y)) a = gadd(a, gatanh(gel(y,2),prec));
     749          14 :       return gerepileupto(av, a);
     750             :   }
     751          21 :   return trans_eval("atanh",gatanh,x,prec);
     752             : }
     753             : /********************************************************************/
     754             : /**                                                                **/
     755             : /**               CACHE BERNOULLI NUMBERS B_2k                     **/
     756             : /**                                                                **/
     757             : /********************************************************************/
     758             : static GEN
     759     1765459 : bern(GEN B, long pr)
     760             : {
     761     1765459 :   if (typ(B) != t_REAL) return fractor(B, pr);
     762      906298 :   if (realprec(B) < pr) return rtor(B,pr);
     763           0 :   return B;
     764             : }
     765             : static const long BERN_MINNB = 5;
     766             : /* need B[2..2*nb] at least prec accuracy. If prec = 0, compute exactly */
     767             : void
     768      187080 : mpbern(long nb, long prec)
     769             : {
     770      187080 :   const pari_sp av = avma;
     771      187080 :   long n, pr, n_is_small = 1, lbern = 0;
     772             :   GEN B;
     773             :   pari_timer T;
     774             : 
     775             :   /* pr = accuracy for computation, prec = required accuracy for result */
     776      187080 :   if (prec)
     777             :   {
     778      185827 :     pr = prec;
     779      185827 :     incrprec(pr);
     780             :   }
     781             :   else
     782        1253 :     pr = prec = LONG_MAX; /* oo */
     783      187080 :   if (nb < BERN_MINNB) nb = BERN_MINNB;
     784      187080 :   if (bernzone)
     785             :   { /* don't recompute known Bernoulli */
     786             :     long i, min, max;
     787      186899 :     lbern = lg(bernzone);
     788      186899 :     if (lbern-1 < nb)
     789         538 :     { min = lbern-1; max = nb; }
     790             :     else
     791      186361 :     { min = nb; max = lbern-1; }
     792             :     /* skip B_2, ..., B_{2*MINNB}, always included as t_FRAC */
     793     7850231 :     for (n = BERN_MINNB+1; n <= min; n++)
     794             :     {
     795     7663522 :       GEN c = gel(bernzone,n);
     796             :       /* also stop if prec = 0 (compute exactly) */
     797     7663522 :       if (typ(c) == t_REAL && realprec(c) < prec) break;
     798             :     }
     799             :     /* B[1..n-1] are OK */
     800      373229 :     if (n > nb) return;
     801         569 :     B = cgetg_block(max+1, t_VEC);
     802         569 :     for (i = 1; i < n; i++) gel(B,i) = gel(bernzone,i);
     803             :     /* keep B[nb+1..max] */
     804         569 :     for (i = nb+1; i <= max; i++) gel(B,i) = gel(bernzone,i);
     805             :   }
     806             :   else
     807             :   {
     808         181 :     B = cgetg_block(nb+1, t_VEC);
     809         181 :     gel(B,1) = gclone(mkfrac(gen_1, utoipos(6)));
     810         181 :     gel(B,2) = gclone(mkfrac(gen_m1, utoipos(30)));
     811         181 :     gel(B,3) = gclone(mkfrac(gen_1, utoipos(42)));
     812         181 :     gel(B,4) = gel(B,2);
     813         181 :     gel(B,5) = gclone(mkfrac(utoipos(5), utoipos(66)));
     814         181 :     n = BERN_MINNB+1;
     815             :   }
     816         750 :   avma = av;
     817         750 :   if (DEBUGLEVEL) {
     818           0 :     err_printf("caching Bernoulli numbers 2 to 2*%ld, prec = %ld\n",
     819             :                nb, prec == LONG_MAX? 0: prec);
     820           0 :     timer_start(&T);
     821             :   }
     822             : 
     823             :   /* B_{2n} = (2n-1) / (4n+2) -
     824             :    * sum_{a = 1}^{n-1} (2n)...(2n+2-2a) / (2...(2a-1)2a) B_{2a} */
     825         750 :   n_is_small = 1;
     826       21144 :   for (; n <= nb; n++, avma = av)
     827             :   { /* compute and store B[n] = B_{2n} */
     828             :     GEN S;
     829       20394 :     if (n < lbern)
     830             :     {
     831        8410 :       GEN b = gel(bernzone,n);
     832        8410 :       if (typ(b)!=t_REAL || realprec(b)>=prec) { gel(B,n) = b; continue; }
     833             :     }
     834             :     /* Not cached, must compute */
     835             :     /* huge accuracy ? May as well compute exactly */
     836       24797 :     if (n_is_small && (prec == LONG_MAX ||
     837        4424 :                        2*n * log((double)2*n) < prec2nbits_mul(prec, LOG2)))
     838        3777 :       S = bernfrac_using_zeta(2*n);
     839             :     else
     840             :     {
     841             : #ifdef LONG_IS_64BIT
     842       13950 :       const ulong mul_overflow = 3037000500UL;
     843             : #else
     844        2646 :       const ulong mul_overflow = 46341UL;
     845             : #endif
     846       16596 :       ulong u = 8, v = 5, a = n-1, b = 2*n-3;
     847       16596 :       n_is_small = 0;
     848       16596 :       S = bern(gel(B,a), pr); /* B_2a */
     849             :       for (;;)
     850             :       { /* b = 2a-1, u = 2v-2, 2a + v = 2n+3 */
     851     1765459 :         if (a == 1) { S = mulri(S, muluu(u,v)); break; } /* a=b=1, v=2n+1, u=4n */
     852             :         /* beware overflow */
     853     1748863 :         S = (v <= mul_overflow)? mulru(S, u*v): mulri(S, muluu(u,v));
     854     1748863 :         S = (a <= mul_overflow)? divru(S, a*b): divri(S, muluu(a,b));
     855     1748863 :         u += 4; v += 2; a--; b -= 2;
     856     1748863 :         S = addrr(bern(gel(B,a), pr), S);
     857     1748863 :         if ((a & 127) == 0) S = gerepileuptoleaf(av, S);
     858     1748863 :       }
     859       16596 :       S = divru(subsr(2*n, S), 2*n+1);
     860       16596 :       shiftr_inplace(S, -2*n);
     861       16596 :       if (realprec(S) != prec) S = rtor(S, prec);
     862             :     }
     863       20373 :     gel(B,n) = gclone(S); /* S = B_2n */
     864             :   }
     865         750 :   if (DEBUGLEVEL) timer_printf(&T, "Bernoulli");
     866         750 :   swap(B, bernzone);
     867         750 :   if (B)
     868             :   { /* kill old non-reused values */
     869       44224 :     for (n = lbern-1; n; n--)
     870             :     {
     871       43655 :       if (gel(B,n) != gel(bernzone,n)) gunclone(gel(B,n));
     872             :     }
     873         569 :     killblock(B);
     874             :   }
     875         750 :   avma = av;
     876             : }
     877             : 
     878             : GEN
     879       11900 : bernfrac(long n)
     880             : {
     881             :   long k;
     882       11900 :   if (n < 0) pari_err_DOMAIN("bernfrac", "index", "<", gen_0, stoi(n));
     883       11893 :   if (n == 0) return gen_1;
     884       10640 :   if (n == 1) return mkfrac(gen_m1,gen_2);
     885        9394 :   if (odd(n)) return gen_0;
     886        8204 :   k = n >> 1;
     887        8204 :   if (!bernzone && k <= BERN_MINNB) mpbern(BERN_MINNB, 0);
     888        8204 :   if (bernzone && k < lg(bernzone))
     889             :   {
     890        6741 :     GEN B = gel(bernzone, k), C;
     891        6741 :     if (typ(B) != t_REAL) return B;
     892        1154 :     C = bernfrac_using_zeta(n);
     893        1154 :     gel(bernzone, k) = gclone(C);
     894        1154 :     gunclone(B); return C;
     895             :   }
     896        1463 :   return bernfrac_using_zeta(n);
     897             : }
     898             : 
     899             : /* mpbern as exact fractions */
     900             : static GEN
     901         168 : bernvec_old(long nb)
     902             : {
     903             :   long n, i;
     904             :   GEN y;
     905             : 
     906         168 :   if (nb < 0) return cgetg(1, t_VEC);
     907          24 :   if (nb > 46340 && BITS_IN_LONG == 32) pari_err_IMPL( "bernvec for n > 46340");
     908             : 
     909         168 :   y = cgetg(nb+2, t_VEC); gel(y,1) = gen_1;
     910         679 :   for (n = 1; n <= nb; n++)
     911             :   { /* compute y[n+1] = B_{2n} */
     912         511 :     pari_sp av = avma;
     913         511 :     GEN b = gmul2n(utoineg(2*n - 1), -1); /* 1 + (2n+1)B_1 = -(2n-1) /2 */
     914         511 :     GEN c = gen_1;
     915         511 :     ulong u1 = 2*n + 1, u2 = n, d1 = 1, d2 = 1;
     916             : 
     917        1113 :     for (i = 1; i < n; i++)
     918             :     {
     919         602 :       c = diviiexact(muliu(c, u1*u2), utoipos(d1*d2));/*= binomial(2n+1, 2*i) */
     920         602 :       b = gadd(b, gmul(c, gel(y,i+1)));
     921         602 :       u1 -= 2; u2--; d1++; d2 += 2;
     922             :     }
     923         511 :     gel(y,n+1) = gerepileupto(av, gdivgs(b, -(1+2*n)));
     924             :   }
     925         168 :   return y;
     926             : }
     927             : GEN
     928         175 : bernvec(long nb)
     929             : {
     930         175 :   long i, l = nb+2;
     931         175 :   GEN y = cgetg(l, t_VEC);
     932         175 :   if (nb < 20) return bernvec_old(nb);
     933           7 :   for (i = 1; i < l; i++) gel(y,i) = bernfrac((i-1) << 1);
     934           7 :   return y;
     935             : }
     936             : 
     937             : /* x := pol_x(v); B_k(x) = \sum_{i=0}^k binomial(k, i) B_i x^{k-i} */
     938             : static GEN
     939        1246 : bernpol_i(long k, long v)
     940             : {
     941             :   GEN B, C;
     942             :   long i;
     943        1246 :   if (v < 0) v = 0;
     944        1246 :   if (k < 0) pari_err_DOMAIN("bernpol", "index", "<", gen_0, stoi(k));
     945        1239 :   mpbern(k >> 1, 0); /* cache B_2, ..., B_2[k/2] */
     946        1239 :   C = vecbinomial(k);
     947        1239 :   B = cgetg(k + 3, t_POL);
     948        1239 :   for (i = 0; i <= k; ++i) gel(B, k-i+2) = gmul(gel(C,i+1), bernfrac(i));
     949        1239 :   B[1] = evalsigne(1) | evalvarn(v);
     950        1239 :   return B;
     951             : }
     952             : GEN
     953        1218 : bernpol(long k, long v) {
     954        1218 :   pari_sp av = avma;
     955        1218 :   return gerepileupto(av, bernpol_i(k, v));
     956             : }
     957             : /* x := pol_x(v); return 1^e + ... + x^e = x^e + (B_{e+1}(x) - B_{e+1})/(e+1) */
     958             : static GEN
     959          42 : faulhaber(long e, long v)
     960             : {
     961             :   GEN B;
     962          42 :   if (e == 0) return pol_x(v);
     963          28 :   B = RgX_integ(bernpol_i(e, v)); /* (B_{e+1}(x) - B_{e+1}) / (e+1) */
     964          28 :   gel(B,e+2) = gaddgs(gel(B,e+2), 1); /* add x^e, in place */
     965          28 :   return B;
     966             : }
     967             : /* sum_v T(v), T a polynomial expression in v */
     968             : GEN
     969          49 : sumformal(GEN T, long v)
     970             : {
     971          49 :   pari_sp av = avma, av2;
     972             :   long i, t, d;
     973             :   GEN R;
     974             : 
     975          49 :   T = simplify_shallow(T);
     976          49 :   t = typ(T);
     977          49 :   if (is_scalar_t(t))
     978          14 :     return gerepileupto(av, monomialcopy(T, 1, v < 0? 0: v));
     979          35 :   if (t != t_POL) pari_err_TYPE("sumformal [not a t_POL]", T);
     980          28 :   if (v < 0) v = varn(T);
     981          28 :   av2 = avma;
     982          28 :   R = gen_0;
     983          28 :   d = poldegree(T,v);
     984          91 :   for (i = d; i >= 0; i--)
     985             :   {
     986          63 :     GEN c = polcoeff0(T, i, v);
     987          63 :     if (gequal0(c)) continue;
     988          42 :     R = gadd(R, gmul(c, faulhaber(i, v)));
     989          42 :     if (gc_needed(av2,3))
     990             :     {
     991           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"sumformal, i = %ld/%ld", i,d);
     992           0 :       R = gerepileupto(av2, R);
     993             :     }
     994             :   }
     995          28 :   return gerepileupto(av, R);
     996             : }
     997             : 
     998             : /********************************************************************/
     999             : /**                                                                **/
    1000             : /**                         EULER'S GAMMA                          **/
    1001             : /**                                                                **/
    1002             : /********************************************************************/
    1003             : 
    1004             : /* x / (i*(i+1)) */
    1005             : GEN
    1006    29118616 : divrunu(GEN x, ulong i)
    1007             : {
    1008    29118616 :   if (i <= LOWMASK) /* i(i+1) < 2^BITS_IN_LONG*/
    1009    29118616 :     return divru(x, i*(i+1));
    1010             :   else
    1011           0 :     return divru(divru(x, i), i+1);
    1012             : }
    1013             : /* x / (i*(i+1)) */
    1014             : GEN
    1015      666137 : divgunu(GEN x, ulong i)
    1016             : {
    1017             : #ifdef LONG_IS_64BIT
    1018      571368 :   if (i < 3037000500L) /* i(i+1) < 2^63 */
    1019             : #else
    1020       94769 :   if (i < 46341L) /* i(i+1) < 2^31 */
    1021             : #endif
    1022      666137 :     return gdivgs(x, i*(i+1));
    1023             :   else
    1024           0 :     return gdivgs(gdivgs(x, i), i+1);
    1025             : }
    1026             : 
    1027             : /* arg(s+it) */
    1028             : double
    1029      181445 : darg(double s, double t)
    1030             : {
    1031             :   double x;
    1032      181445 :   if (!t) return (s>0)? 0.: M_PI;
    1033      171675 :   if (!s) return (t>0)? M_PI/2: -M_PI/2;
    1034      171668 :   x = atan(t/s);
    1035      171668 :   return (s>0)? x
    1036      171668 :               : ((t>0)? x+M_PI : x-M_PI);
    1037             : }
    1038             : 
    1039             : void
    1040      181445 : dcxlog(double s, double t, double *a, double *b)
    1041             : {
    1042      181445 :   *a = log(s*s + t*t) / 2; /* log |s| = Re(log(s)) */
    1043      181445 :   *b = darg(s,t);          /* Im(log(s)) */
    1044      181445 : }
    1045             : 
    1046             : double
    1047       12838 : dabs(double s, double t) { return sqrt( s*s + t*t ); }
    1048             : double
    1049          35 : dnorm(double s, double t) { return s*s + t*t; }
    1050             : 
    1051             : #if 0
    1052             : /* x, z t_REAL. Compute unique x in ]-z,z] congruent to x mod 2z */
    1053             : static GEN
    1054             : red_mod_2z(GEN x, GEN z)
    1055             : {
    1056             :   GEN Z = gmul2n(z, 1), d = subrr(z, x);
    1057             :   /* require little accuracy */
    1058             :   if (!signe(d)) return x;
    1059             :   setprec(d, nbits2prec(expo(d) - expo(Z)));
    1060             :   return addrr(mulir(floorr(divrr(d, Z)), Z), x);
    1061             : }
    1062             : #endif
    1063             : 
    1064             : /* lngamma(1+z) = -Euler*z + sum_{i > 1} zeta(i)/i (-z)^i
    1065             :  * at relative precision prec, |z| < 1 is small */
    1066             : static GEN
    1067        2030 : lngamma1(GEN z, long prec)
    1068             : { /* sum_{i > l} |z|^(i-1) = |z|^l / (1-|z|) < 2^-B
    1069             :    * for l > (B+1) / |log2(|z|)| */
    1070        2030 :   long i, l = ceil((bit_accuracy(prec) + 1) / - dbllog2(z));
    1071        2030 :   GEN zet, me = mpeuler(prec), s = gen_0;
    1072        2030 :   setsigne(me, -1); /* -Euler */
    1073        2030 :   if (l <= 1) return gmul(me, z);
    1074        1855 :   zet = veczeta(gen_1, gen_2, l-1, prec); /* z[i] = zeta(i+1) */
    1075       15370 :   for (i = l; i > 1; i--)
    1076             :   {
    1077       13515 :     GEN c = divru(gel(zet,i-1), i);
    1078       13515 :     if (odd(i)) setsigne(c, -1);
    1079       13515 :     s = gadd(gmul(s,z), c);
    1080             :   }
    1081        1855 :   return gmul(z, gadd(gmul(s,z), me));
    1082             : }
    1083             : 
    1084             : static GEN
    1085      181291 : cxgamma(GEN s0, int dolog, long prec)
    1086             : {
    1087             :   GEN s, u, a, y, res, tes, sig, tau, invn2, p1, nnx, pi, pi2, sqrtpi2;
    1088             :   long i, lim, nn, esig, et;
    1089             :   pari_sp av, av2;
    1090      181291 :   int funeq = 0;
    1091             :   pari_timer T;
    1092             : 
    1093      181291 :   if (DEBUGLEVEL>5) timer_start(&T);
    1094      181291 :   s = trans_fix_arg(&prec,&s0,&sig,&tau,&av,&res);
    1095             : 
    1096      181291 :   esig = expo(sig);
    1097      181291 :   et = signe(tau)? expo(tau): 0;
    1098      181291 :   if ((signe(sig) <= 0 || esig < -1) && et <= 16)
    1099             :   { /* s <--> 1-s */
    1100        3045 :     funeq = 1; s = gsubsg(1, s); sig = real_i(s);
    1101             :   }
    1102             : 
    1103             :   /* find "optimal" parameters [lim, nn] */
    1104      181291 :   if (esig > 300 || et > 300)
    1105          35 :   { /* |s| is HUGE ! Play safe and avoid inf / NaN */
    1106             :     GEN S, iS, l2, la, u;
    1107             :     double logla, l;
    1108             : 
    1109          35 :     S = gprec_w(s,LOWDEFAULTPREC);
    1110             :     /* l2 ~ |lngamma(s))|^2 */
    1111          35 :     l2 = gnorm(gmul(S, glog(S, LOWDEFAULTPREC)));
    1112          35 :     l = (prec2nbits_mul(prec, LOG2) - rtodbl(glog(l2,LOWDEFAULTPREC))/2) / 2.;
    1113          35 :     if (l < 0) l = 0.;
    1114             : 
    1115          35 :     iS = imag_i(S);
    1116          35 :     if (et > 0 && l > 0)
    1117          21 :     {
    1118          21 :       GEN t = gmul(iS, dbltor(M_PI / l)), logt = glog(t,LOWDEFAULTPREC);
    1119          21 :       la = gmul(t, logt);
    1120          21 :       if      (gcmpgs(la, 3) < 0)   { logla = log(3.); la = stoi(3); }
    1121          14 :       else if (gcmpgs(la, 150) > 0) { logla = rtodbl(logt); la = t; }
    1122           7 :       else logla = rtodbl(mplog(la));
    1123             :     }
    1124             :     else
    1125             :     {
    1126          14 :       logla = log(3.); la = stoi(3);
    1127             :     }
    1128          35 :     lim = (long)ceil(l / (1.+ logla));
    1129          35 :     if (lim == 0) lim = 1;
    1130             : 
    1131          35 :     u = gmul(la, dbltor((lim-0.5)/M_PI));
    1132          35 :     l2 = gsub(gsqr(u), gsqr(iS));
    1133          35 :     if (signe(l2) > 0)
    1134             :     {
    1135          14 :       l2 = gsub(gsqrt(l2,3), sig);
    1136          14 :       if (signe(l2) > 0) nn = itos( gceil(l2) ); else nn = 1;
    1137             :     }
    1138             :     else
    1139          21 :       nn = 1;
    1140             :   }
    1141             :   else
    1142             :   { /* |s| is moderate. Use floats  */
    1143      181256 :     double ssig = rtodbl(sig);
    1144      181256 :     double st = typ(s) == t_REAL? 0.0: rtodbl(imag_i(s));
    1145             :     double la, l,l2,u,v, rlogs, ilogs;
    1146             : 
    1147      181256 :     if (fabs(ssig-1) + fabs(st) < 1e-16)
    1148             :     { /* s ~ 1: loggamma(1+u) ~ - Euler * u, cancellation */
    1149        2436 :       if (funeq) /* s0 ~ 0: use lngamma(s0)+log(s0) = lngamma(s0+1) */
    1150             :       {
    1151           0 :         if (dolog)
    1152           0 :           y = gsub(lngamma1(s0,prec), glog(s0,prec));
    1153             :         else
    1154           0 :           y = gdiv(gexp(lngamma1(s0,prec), prec), s0);
    1155             :       }
    1156             :       else
    1157             :       {
    1158        4872 :         if (isint1(s0)) { avma = av; return dolog? real_0(prec): real_1(prec); }
    1159         805 :         y = lngamma1(gsubgs(s0,1),prec);
    1160         805 :         if (!dolog) y = gexp(y,prec);
    1161             :       }
    1162         805 :       avma = av; return affc_fixlg(y, res);
    1163             :     }
    1164      178820 :     dcxlog(ssig,st, &rlogs,&ilogs);
    1165             :     /* Re (s - 1/2) log(s) */
    1166      178820 :     u = (ssig - 0.5)*rlogs - st * ilogs;
    1167             :     /* Im (s - 1/2) log(s) */
    1168      178820 :     v = (ssig - 0.5)*ilogs + st * rlogs;
    1169             :     /* l2 = | (s - 1/2) log(s) - s + log(2Pi)/2 |^2 ~ |lngamma(s))|^2 */
    1170      178820 :     u = u - ssig + log(2.*M_PI)/2;
    1171      178820 :     v = v - st;
    1172      178820 :     l2 = u*u + v*v;
    1173      178820 :     if (l2 < 0.000001) l2 = 0.000001;
    1174      178820 :     l = (prec2nbits_mul(prec, LOG2) - log(l2)/2) / 2.;
    1175      178820 :     if (l < 0) l = 0.;
    1176             : 
    1177      178820 :     la = 3.; /* FIXME: heuristic... */
    1178      178820 :     if (st > 1 && l > 0)
    1179             :     {
    1180       63378 :       double t = st * M_PI / l;
    1181       63378 :       la = t * log(t);
    1182       63378 :       if (la < 3) la = 3.;
    1183       63378 :       if (la > 150) la = t;
    1184             :     }
    1185      178820 :     lim = (long)ceil(l / (1.+ log(la)));
    1186      178820 :     if (lim == 0) lim = 1;
    1187             : 
    1188      178820 :     u = (lim-0.5) * la / M_PI;
    1189      178820 :     l2 = u*u - st*st;
    1190      178820 :     if (l2 > 0)
    1191             :     {
    1192      164727 :       double t = ceil(sqrt(l2) - ssig);
    1193      164727 :       nn = (t < 1)? 1: (long)t;
    1194      164727 :       if (nn < 1) nn = 1;
    1195             :     }
    1196             :     else
    1197       14093 :       nn = 1;
    1198      178820 :     if (DEBUGLEVEL>5) err_printf("lim, nn: [%ld, %ld], la = %lf\n",lim,nn,la);
    1199             :   }
    1200      178855 :   incrprec(prec);
    1201             : 
    1202      178855 :   av2 = avma;
    1203      178855 :   y = s;
    1204      178855 :   if (typ(s0) == t_INT)
    1205             :   {
    1206        1748 :     if (signe(s0) <= 0)
    1207           0 :       pari_err_DOMAIN("gamma","argument", "=",
    1208             :                        strtoGENstr("non-positive integer"), s0);
    1209        1748 :     if (is_bigint(s0)) {
    1210           7 :       for (i=1; i < nn; i++)
    1211             :       {
    1212           0 :         y = mulri(y, addiu(s0, i));
    1213           0 :         if (gc_needed(av2,3))
    1214             :         {
    1215           0 :           if(DEBUGMEM>1) pari_warn(warnmem,"gamma");
    1216           0 :           y = gerepileuptoleaf(av2, y);
    1217             :         }
    1218             :       }
    1219             :     } else {
    1220        1741 :       ulong ss = itou(s0);
    1221       16470 :       for (i=1; i < nn; i++)
    1222             :       {
    1223       14729 :         y = mulru(y, ss + i);
    1224       14729 :         if (gc_needed(av2,3))
    1225             :         {
    1226           0 :           if(DEBUGMEM>1) pari_warn(warnmem,"gamma");
    1227           0 :           y = gerepileuptoleaf(av2, y);
    1228             :         }
    1229             :       }
    1230             :     }
    1231        1748 :     if (dolog) y = logr_abs(y);
    1232             :   }
    1233             :   else
    1234             :   { /* Compute lngamma mod 2 I Pi */
    1235      177107 :     GEN sq = gsqr(s);
    1236      177107 :     pari_sp av3 = avma;
    1237     3306260 :     for (i = 1; i < nn - 1; i += 2)
    1238             :     {
    1239     3129153 :       y = gmul(y, gaddsg(i*(i + 1), gadd(gmulsg(2*i + 1, s), sq)));
    1240     3129153 :       if (gc_needed(av2,3))
    1241             :       {
    1242           0 :         if(DEBUGMEM>1) pari_warn(warnmem,"gamma");
    1243           0 :         y = gerepileupto(av3, y);
    1244             :       }
    1245             :     }
    1246      177107 :     if (!odd(nn)) y = gmul(y, gaddsg(nn - 1, s));
    1247      177107 :     if (dolog)
    1248             :     {
    1249        6321 :       if (typ(s) == t_REAL) y = logr_abs(y);
    1250             :       else
    1251             :       { /* fix imaginary part */
    1252         273 :         long prec0 = LOWDEFAULTPREC;
    1253         273 :         GEN s0 = gprec_w(s, prec0), y0 = s0, k;
    1254         273 :         y0 = garg(y0, prec0); /* Im log(s) at low accuracy */
    1255         273 :         for (i=1; i < nn; i++) y0 = gadd(y0, garg(gaddgs(s0,i), prec0));
    1256         273 :         y = glog(y, prec);
    1257         273 :         k = ground( gdiv(gsub(y0, imag_i(y)), Pi2n(1,prec0)) );
    1258         273 :         if (signe(k)) y = gadd(y, mulcxI(mulir(k, Pi2n(1, prec))));
    1259             :       }
    1260             :     }
    1261             :   }
    1262      178855 :   if (DEBUGLEVEL>5) timer_printf(&T,"product from 0 to N-1");
    1263             : 
    1264      178855 :   nnx = gaddgs(s, nn);
    1265      178855 :   a = ginv(nnx); invn2 = gsqr(a);
    1266      178855 :   av2 = avma;
    1267      178855 :   mpbern(lim,prec);
    1268      178855 :   tes = divrunu(bernreal(2*lim,prec), 2*lim-1); /* B2l / (2l-1) 2l*/
    1269      178855 :   if (DEBUGLEVEL>5) timer_printf(&T,"Bernoullis");
    1270     7934203 :   for (i = 2*lim-2; i > 1; i -= 2)
    1271             :   {
    1272     7755348 :     u = divrunu(bernreal(i,prec), i-1); /* Bi / i(i-1) */
    1273     7755348 :     tes = gadd(u, gmul(invn2,tes));
    1274     7755348 :     if (gc_needed(av2,3))
    1275             :     {
    1276           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"gamma");
    1277           0 :       tes = gerepileupto(av2, tes);
    1278             :     }
    1279             :   }
    1280      178855 :   if (DEBUGLEVEL>5) timer_printf(&T,"Bernoulli sum");
    1281             : 
    1282      178855 :   p1 = gsub(gmul(gsub(nnx, ghalf), glog(nnx,prec)), nnx);
    1283      178855 :   p1 = gadd(p1, gmul(tes, a));
    1284             : 
    1285      178855 :   pi = mppi(prec); pi2 = shiftr(pi, 1); sqrtpi2 = sqrtr(pi2);
    1286             : 
    1287      178855 :   if (dolog)
    1288             :   {
    1289        6354 :     if (funeq)
    1290             :     { /* recall that s = 1 - s0 */
    1291          14 :       GEN T = shiftr(sqrtpi2,-1); /* sqrt(2Pi)/2 */
    1292          14 :       if (typ(s) != t_REAL)
    1293             :       {
    1294             :         /* We compute log(sin(Pi s0)) so that it has branch cuts along
    1295             :         * (-oo, 0] and [1, oo).  To do this in a numerically stable way
    1296             :         * we must compute the log first then mangle its imaginary part.
    1297             :         * The rounding operation below is stable because we're rounding
    1298             :         * a number which is already within 1/4 of an integer. */
    1299             : 
    1300             :         /* z = log( sin(Pi s0) / (sqrt(2Pi)/2) ) */
    1301           7 :         GEN z = glog(gdiv(gsin(gmul(pi,s0),prec), T), prec);
    1302             :         /* b = (2 Re(s) - 1) / 4 */
    1303           7 :         GEN b = shiftr(subrs(shiftr(sig, 1), 1), -2);
    1304           7 :         y = gsub(y, z);
    1305           7 :         if (gsigne(imag_i(s)) > 0) togglesign(b);
    1306             :         /* z = 2Pi round( Im(z)/2Pi - b ) */
    1307           7 :         z = gmul(roundr(gsub(gdiv(imag_i(z), pi2), b)), pi2);
    1308           7 :         if (signe(z)) { /* y += I*z, z a t_REAL */
    1309           0 :           if (typ(y) == t_COMPLEX)
    1310           0 :             gel(y,2) = gadd(gel(y,2), z);
    1311             :           else
    1312           0 :             y = mkcomplex(y, z);
    1313             :         }
    1314             :       }
    1315             :       else
    1316             :       { /* s0 < 0, formula simplifies: imag(lngamma(s0)) = - Pi * floor(s0) */
    1317           7 :         GEN z = logr_abs(divrr(mpsin(gmul(pi,s0)), T));
    1318           7 :         y = gsub(y, z);
    1319           7 :         y = mkcomplex(y, mulri(pi, gfloor(s0)));
    1320             :       }
    1321          14 :       p1 = gneg(p1);
    1322             :     }
    1323             :     else /* y --> sqrt(2Pi) / y */
    1324        6340 :       y = gsub(logr_abs(sqrtpi2), y);
    1325        6354 :     y = gadd(p1, y);
    1326             :   }
    1327             :   else
    1328             :   {
    1329      172501 :     if (funeq)
    1330             :     { /* y --> y Pi/(sin(Pi s) * sqrt(2Pi)) = y sqrt(Pi/2)/sin(Pi s) */
    1331        3031 :       y = gdiv(gmul(shiftr(sqrtpi2,-1),y), gsin(gmul(pi,s0), prec));
    1332             :       /* don't use s above: sin(pi s0) = sin(pi s) and the former is
    1333             :        * more accurate, esp. if s0 ~ 0 */
    1334        3031 :       p1 = gneg(p1);
    1335             :     }
    1336             :     else /* y --> sqrt(2Pi) / y */
    1337      169470 :       y = gdiv(sqrtpi2, y);
    1338      172501 :     y = gmul(gexp(p1, prec), y);
    1339             :   }
    1340      178855 :   avma = av; return affc_fixlg(y, res);
    1341             : }
    1342             : 
    1343             : /* Gamma((m+1) / 2) */
    1344             : static GEN
    1345       10913 : gammahs(long m, long prec)
    1346             : {
    1347       10913 :   GEN y = cgetr(prec), z;
    1348       10913 :   pari_sp av = avma;
    1349       10913 :   long ma = labs(m);
    1350             : 
    1351       10913 :   if (ma > 200 + 50*(prec-2)) /* heuristic */
    1352             :   {
    1353           7 :     z = stor(m + 1, prec); shiftr_inplace(z, -1);
    1354           7 :     affrr(cxgamma(z,0,prec), y);
    1355           7 :     avma = av; return y;
    1356             :   }
    1357       10906 :   z = sqrtr( mppi(prec) );
    1358       10906 :   if (m)
    1359             :   {
    1360        3087 :     GEN p1 = mulu_interval(ma/2 + 1, ma);
    1361        3087 :     long v = vali(p1);
    1362        3087 :     p1 = shifti(p1, -v); v -= ma;
    1363        3087 :     if (m >= 0) z = mulri(z,p1);
    1364             :     else
    1365             :     {
    1366          63 :       z = divri(z,p1); v = -v;
    1367          63 :       if ((m&3) == 2) setsigne(z,-1);
    1368             :     }
    1369        3087 :     shiftr_inplace(z, v);
    1370             :   }
    1371       10906 :   affrr(z, y); avma = av; return y;
    1372             : }
    1373             : GEN
    1374          21 : ggammah(GEN x, long prec)
    1375             : {
    1376          21 :   switch(typ(x))
    1377             :   {
    1378             :     case t_INT:
    1379             :     {
    1380          14 :       long k = itos(x);
    1381          14 :       if (labs(k) > 962353) pari_err_OVERFLOW("gammah");
    1382          14 :       return gammahs(k<<1, prec);
    1383             :     }
    1384             :     case t_REAL: case t_COMPLEX: case t_PADIC: case t_SER: {
    1385           7 :       pari_sp av = avma;
    1386           7 :       return gerepileupto(av, ggamma(gadd(x,ghalf), prec));
    1387             :     }
    1388             :   }
    1389           0 :   return trans_eval("gammah",ggammah,x,prec);
    1390             : }
    1391             : 
    1392             : /* find n such that n+v_p(n!)>=k p^2/(p-1)^2 */
    1393             : static long
    1394          28 : nboft(long k, long p)
    1395             : {
    1396          28 :   pari_sp av = avma;
    1397             :   long s, n;
    1398             : 
    1399          28 :   if (k <= 0) return 0;
    1400          28 :   k = itou( gceil(gdiv(mului(k, sqru(p)), sqru(p-1))) );
    1401          28 :   avma = av;
    1402          28 :   for (s=0, n=0; n+s < k; n++, s += u_lval(n, p));
    1403          28 :   return n;
    1404             : }
    1405             : 
    1406             : /* Using Dwork's expansion, compute \Gamma(px+1)=-\Gamma(px) with x a unit.
    1407             :  * See p-Adic Gamma Functions and Dwork Cohomology, Maurizio Boyarsky
    1408             :  * Transactions of the AMS, Vol. 257, No. 2. (Feb., 1980), pp. 359-369.
    1409             :  * Inspired by a GP script by Fernando Rodriguez-Villegas */
    1410             : static GEN
    1411          28 : gadw(GEN x, long p)
    1412             : {
    1413          28 :   pari_sp ltop = avma;
    1414          28 :   GEN s, t, u = cgetg(p+1, t_VEC);
    1415          28 :   long j, k, kp, n = nboft(precp(x)+valp(x)+1, p);
    1416             : 
    1417          28 :   t = s = gaddsg(1, zeropadic(gel(x,2), n));
    1418          28 :   gel(u, 1) = s;
    1419          28 :   gel(u, 2) = s;
    1420          42 :   for (j = 2; j < p; ++j)
    1421          14 :     gel(u, j+1) = gdivgs(gel(u, j), j);
    1422         378 :   for (k = 1, kp = p; k < n; ++k, kp += p) /* kp = k*p */
    1423             :   {
    1424             :     GEN c;
    1425         350 :     gel(u, 1) = gdivgs(gadd(gel(u, 1), gel(u, p)), kp);
    1426         812 :     for (j = 1; j < p; ++j)
    1427         462 :       gel(u, j+1) = gdivgs(gadd(gel(u, j), gel(u, j+1)), kp + j);
    1428             : 
    1429         350 :     t = gmul(t, gaddgs(x, k-1));
    1430         350 :     c = leafcopy(gel(u,1)); setvalp(c, valp(c) + k); /* c = u[1] * p^k */
    1431         350 :     s = gadd(s, gmul(c, t));
    1432         350 :     if ((k&0xFL)==0) gerepileall(ltop, 3, &u,&s,&t);
    1433             :   }
    1434          28 :   return gneg(s);
    1435             : }
    1436             : 
    1437             : /*Use Dwork expansion*/
    1438             : /*This is a O(p*e*log(pe)) algorithm, should be used when p small
    1439             :  * If p==2 this is a O(pe) algorithm. */
    1440             : static GEN
    1441          28 : Qp_gamma_Dwork(GEN x, long p)
    1442             : {
    1443          28 :   pari_sp ltop = avma;
    1444          28 :   long k = padic_to_Fl(x, p);
    1445             :   GEN p1;
    1446             :   long j;
    1447          28 :   long px = precp(x);
    1448          28 :   if (p==2 && px)
    1449             :   {
    1450          14 :     x = shallowcopy(x);
    1451          14 :     setprecp(x, px+1);
    1452          14 :     gel(x,3) = shifti(gel(x,3),1);
    1453             :   }
    1454          28 :   if (k)
    1455             :   {
    1456          21 :     GEN x_k = gsubgs(x,k);
    1457          21 :     x = gdivgs(x_k, p);
    1458          21 :     p1 = gadw(x, p); if (!odd(k)) p1 = gneg(p1);
    1459          21 :     for (j = 1; j < k; ++j) p1 = gmul(p1, gaddgs(x_k, j));
    1460             :   }
    1461             :   else
    1462           7 :     p1 = gneg(gadw(gdivgs(x, p), p));
    1463          28 :   return gerepileupto(ltop, p1);
    1464             : }
    1465             : 
    1466             : /* Compute Qp_gamma using the definition. This is a O(x*M(log(pe))) algorithm.
    1467             :  * This should be used if x is very small. */
    1468             : static GEN
    1469          49 : Qp_gamma_Morita(long n, GEN p, long e)
    1470             : {
    1471          49 :   pari_sp ltop=avma;
    1472          49 :   GEN p2 = gaddsg((n&1)?-1:1, zeropadic(p, e));
    1473             :   long i;
    1474          49 :   long pp=is_bigint(p)? 0: itos(p);
    1475         154 :   for (i = 2; i < n; i++)
    1476         105 :     if (!pp || i%pp)
    1477             :     {
    1478          63 :       p2 = gmulgs(p2, i);
    1479          63 :       if ((i&0xFL) == 0xFL)
    1480           0 :         p2 = gerepileupto(ltop, p2);
    1481             :     }
    1482          49 :   return gerepileupto(ltop, p2);
    1483             : }
    1484             : 
    1485             : /* x\in\N: Gamma(-x)=(-1)^(1+x+x\p)*Gamma(1+x) */
    1486             : static GEN
    1487          28 : Qp_gamma_neg_Morita(long n, GEN p, long e)
    1488             : {
    1489          28 :   GEN g = ginv(Qp_gamma_Morita(n+1, p, e));
    1490          28 :   return ((n^sdivsi(n,p)) & 1)? g: gneg(g);
    1491             : }
    1492             : 
    1493             : /* p-adic Gamma function for x a p-adic integer */
    1494             : /* If n < p*e : use Morita's definition.
    1495             :  * Else : use Dwork's expansion.
    1496             :  * If both n and p are big : itos(p) will fail.
    1497             :  * TODO: handle p=2 better (Qp_gamma_Dwork is slow for p=2). */
    1498             : GEN
    1499          77 : Qp_gamma(GEN x)
    1500             : {
    1501          77 :   GEN n, m, N, p = gel(x,2);
    1502          77 :   long s, e = precp(x);
    1503          77 :   if (absequaliu(p, 2) && e == 2) e = 1;
    1504          77 :   if (valp(x) < 0) pari_err_DOMAIN("gamma","v_p(x)", "<", gen_0, x);
    1505          77 :   n = gtrunc(x);
    1506          77 :   m = gtrunc(gneg(x));
    1507          77 :   N = cmpii(n,m)<=0?n:m;
    1508          77 :   s = itos_or_0(N);
    1509          77 :   if (s && cmpsi(s, muliu(p,e)) < 0) /* s < p*e */
    1510          49 :     return (N == n) ? Qp_gamma_Morita(s,p,e): Qp_gamma_neg_Morita(s,p,e);
    1511          28 :   return Qp_gamma_Dwork(x, itos(p));
    1512             : }
    1513             : 
    1514             : /* gamma(1+x) - 1, |x| < 1 is "small" */
    1515             : GEN
    1516        1211 : ggamma1m1(GEN x, long prec) { return gexpm1(lngamma1(x, prec), prec); }
    1517             : 
    1518             : /* lngamma(y) with 0 constant term, using (lngamma y)' = y' psi(y) */
    1519             : static GEN
    1520       16016 : serlngamma0(GEN y, long prec)
    1521             : {
    1522             :   GEN t;
    1523       16016 :   if (valp(y)) pari_err_DOMAIN("lngamma","valuation", "!=", gen_0, y);
    1524       16009 :   t = derivser(y);
    1525             :   /* don't compute psi if y'=0 */
    1526       16009 :   if (signe(t)) t = gmul(t, gpsi(y,prec));
    1527       16009 :   return integser(t);
    1528             : }
    1529             : 
    1530             : static GEN
    1531       15981 : serlngamma(GEN y, long prec)
    1532             : {
    1533             :   GEN z, y0, Y;
    1534       15981 :   if (lg(y) == 2) pari_err_DOMAIN("gamma", "argument", "=", gen_0,y);
    1535             :   /* exp(lngamma) */
    1536       15974 :   if (valp(y) > 0) return gdiv(gexp(glngamma(gaddgs(y,1),prec),prec),y);
    1537       15757 :   y0 = simplify_shallow(gel(y,2));
    1538       15757 :   z = NULL; Y = y;
    1539       15757 :   if (isint(y0, &y0))
    1540             :   { /* fun eq. avoids log singularity of lngamma at negative ints */
    1541        6300 :     long s = signe(y0);
    1542             :     /* possible if y[2] is an inexact 0 */
    1543        6300 :     if (!s) return gdiv(gexp(glngamma(gaddgs(y,1),prec),prec),y);
    1544        6293 :     if (signe(y0) < 0) { Y = gsubsg(1, y); y0 = subsi(1, y0); }
    1545        6293 :     if (abscmpiu(y0, 50) < 0) z = mpfact(itos(y0)-1); /* more precise */
    1546             :   }
    1547       15750 :   if (!z) z = ggamma(y0,prec);
    1548       15750 :   z = gmul(z, gexp(serlngamma0(Y,prec),prec));
    1549       15750 :   if (Y != y)
    1550             :   {
    1551          77 :     GEN pi = mppi(prec);
    1552          77 :     z = gdiv(mpodd(y0)? pi: negr(pi),
    1553             :              gmul(z, gsin(gmul(pi,serchop0(y)), prec)));
    1554             :   }
    1555       15750 :   return z;
    1556             : }
    1557             : GEN
    1558      206423 : ggamma(GEN x, long prec)
    1559             : {
    1560             :   pari_sp av;
    1561             :   GEN y;
    1562             : 
    1563      206423 :   switch(typ(x))
    1564             :   {
    1565             :     case t_INT:
    1566        5369 :       if (signe(x) <= 0)
    1567           0 :         pari_err_DOMAIN("gamma","argument", "=",
    1568             :                          strtoGENstr("non-positive integer"), x);
    1569        5369 :       if (abscmpiu(x,481177) > 0) pari_err_OVERFLOW("gamma");
    1570        5369 :       return mpfactr(itos(x) - 1, prec);
    1571             : 
    1572             :     case t_REAL: case t_COMPLEX:
    1573      173551 :       return cxgamma(x, 0, prec);
    1574             : 
    1575             :     case t_FRAC:
    1576             :     {
    1577       11452 :       GEN a = gel(x,1), b = gel(x,2), c;
    1578             :       long m;
    1579       11452 :       if (absequaliu(b,2))
    1580             :       {
    1581       10899 :         if (is_bigint(a) || labs(m = itos(a)) > 962354)
    1582             :         {
    1583           0 :           pari_err_OVERFLOW("gamma");
    1584             :           return NULL; /* LCOV_EXCL_LINE */
    1585             :         }
    1586       10899 :         return gammahs(m-1, prec);
    1587             :       }
    1588         553 :       av = avma; c = subii(a,b);
    1589         553 :       if (expi(c) - expi(b) < -50)
    1590             :       { /* x = 1 + c/b is close to 1 */
    1591           7 :         x = mkfrac(c,b);
    1592           7 :         if (lgefint(b) >= prec) x = fractor(x,prec);
    1593           7 :         y = mpexp(lngamma1(x, prec));
    1594             :       }
    1595         546 :       else if (signe(a) < 0 || cmpii(shifti(a,1), b) < 0)
    1596         217 :       { /* gamma will use functional equation x -> z = 1-x = -c/b >= 1/2.
    1597             :          * Gamma(x) = Pi / (sin(Pi z) * Gamma(z)) */
    1598         217 :         GEN z = mkfrac(negi(c), b), q = ground(z), r = gsub(z,q);
    1599         217 :         GEN pi = mppi(prec); /* |r| <= 1/2 */
    1600         217 :         z = fractor(z, prec+EXTRAPRECWORD);
    1601         217 :         y = divrr(pi, mulrr(mpsin(gmul(pi, r)), cxgamma(z, 0, prec)));
    1602         217 :         if (mpodd(q)) togglesign(y);
    1603             :       }
    1604             :       else
    1605             :       {
    1606         329 :         x = fractor(x, prec);
    1607         329 :         y = cxgamma(x, 0, prec);
    1608             :       }
    1609         553 :       return gerepileupto(av, y);
    1610             :     }
    1611             : 
    1612          70 :     case t_PADIC: return Qp_gamma(x);
    1613             :     default:
    1614       15981 :       av = avma; if (!(y = toser_i(x))) break;
    1615       15981 :       return gerepileupto(av, serlngamma(y, prec));
    1616             :   }
    1617           0 :   return trans_eval("gamma",ggamma,x,prec);
    1618             : }
    1619             : 
    1620             : GEN
    1621       19260 : mpfactr(long n, long prec)
    1622             : {
    1623       19260 :   GEN f = cgetr(prec);
    1624       19260 :   pari_sp av = avma;
    1625             : 
    1626       19260 :   if (n+1 > 350 + 70*(prec-2)) /* heuristic */
    1627          21 :     affrr(cxgamma(stor(n+1, prec), 0, prec), f);
    1628             :   else
    1629       19239 :     affir(mpfact(n), f);
    1630       19260 :   avma = av; return f;
    1631             : }
    1632             : 
    1633             : GEN
    1634        7595 : glngamma(GEN x, long prec)
    1635             : {
    1636        7595 :   pari_sp av = avma;
    1637             :   GEN y, y0, t;
    1638             : 
    1639        7595 :   switch(typ(x))
    1640             :   {
    1641             :     case t_INT:
    1642         161 :       if (signe(x) <= 0)
    1643           0 :         pari_err_DOMAIN("lngamma","argument", "=",
    1644             :                          strtoGENstr("non-positive integer"), x);
    1645         161 :       if (abscmpiu(x,200 + 50*(prec-2)) > 0) /* heuristic */
    1646          19 :         return cxgamma(x, 1, prec);
    1647         142 :       return gerepileuptoleaf(av, logr_abs( itor(mpfact(itos(x) - 1), prec) ));
    1648             :     case t_FRAC:
    1649             :     {
    1650         903 :       GEN a = gel(x,1), b = gel(x,2), c = subii(a,b);
    1651         903 :       long e = expi(b) - expi(c);
    1652         903 :       if (e > 50)
    1653             :       {
    1654           7 :         x = mkfrac(c,b);
    1655           7 :         if (lgefint(b) >= prec) x = fractor(x,prec + nbits2nlong(e));
    1656           7 :         y = lngamma1(x, prec);
    1657             :       }
    1658         896 :       else if (signe(a) < 0 || cmpii(shifti(a,1), b) < 0)
    1659          21 :       { /* gamma will use functional equation x -> z = 1-x = -c/b >= 1/2.
    1660             :          * lngamma(x) = log |Pi / (sin(Pi z) * Gamma(z))| + I*Pi * floor(x) */
    1661          21 :         GEN z = mkfrac(negi(c), b), q = ground(z), r = gsub(z,q);
    1662          21 :         GEN pi = mppi(prec); /* |r| <= 1/2 */
    1663          21 :         z = fractor(z, prec+EXTRAPRECWORD);
    1664          21 :         y = subrr(logr_abs(divrr(pi, mpsin(gmul(pi, r)))), cxgamma(z, 1, prec));
    1665          21 :         if (signe(a) < 0) y = gadd(y, mkcomplex(gen_0, mulri(pi, gfloor(x))));
    1666             :       }
    1667             :       else
    1668             :       {
    1669         875 :         x = fractor(x, e > 1? prec+EXTRAPRECWORD: prec);
    1670         875 :         y = cxgamma(x, 1, prec);
    1671             :       }
    1672         903 :       return gerepileupto(av, y);
    1673             :     }
    1674             : 
    1675             :     case t_REAL: case t_COMPLEX:
    1676        6251 :       return cxgamma(x, 1, prec);
    1677             : 
    1678             :     default:
    1679         273 :       if (!(y = toser_i(x))) break;
    1680         273 :       if (lg(y) == 2) pari_err_DOMAIN("lngamma", "argument", "=", gen_0,y);
    1681         266 :       t = serlngamma0(y,prec);
    1682         252 :       y0 = simplify_shallow(gel(y,2));
    1683             :       /* no constant term if y0 = 1 or 2 */
    1684         252 :       if (!isint(y0,&y0) || signe(y0) <= 0 || abscmpiu(y0,2) > 2)
    1685           7 :         t = gadd(t, glngamma(y0,prec));
    1686         252 :       return gerepileupto(av, t);
    1687             : 
    1688           7 :     case t_PADIC: return gerepileupto(av, Qp_log(Qp_gamma(x)));
    1689             :   }
    1690           0 :   return trans_eval("lngamma",glngamma,x,prec);
    1691             : }
    1692             : /********************************************************************/
    1693             : /**                                                                **/
    1694             : /**                  PSI(x) = GAMMA'(x)/GAMMA(x)                   **/
    1695             : /**                                                                **/
    1696             : /********************************************************************/
    1697             : static GEN
    1698          42 : cxpsi(GEN s0, long prec)
    1699             : {
    1700             :   pari_sp av, av2;
    1701             :   GEN sum,z,a,res,tes,in2,sig,tau,s,unr,s2,sq;
    1702             :   long lim,nn,k;
    1703          42 :   const long la = 3;
    1704          42 :   int funeq = 0;
    1705             :   pari_timer T;
    1706             : 
    1707          42 :   if (DEBUGLEVEL>2) timer_start(&T);
    1708          42 :   s = trans_fix_arg(&prec,&s0,&sig,&tau,&av,&res);
    1709          42 :   if (signe(sig) <= 0) { funeq = 1; s = gsub(gen_1, s); sig = real_i(s); }
    1710          42 :   if (typ(s0) == t_INT && signe(s0) <= 0)
    1711           0 :     pari_err_DOMAIN("psi","argument", "=",
    1712             :                     strtoGENstr("non-positive integer"), s0);
    1713             : 
    1714          42 :   if (expo(sig) > 300 || (typ(s) == t_COMPLEX && gexpo(gel(s,2)) > 300))
    1715           7 :   { /* |s| is HUGE. Play safe */
    1716           7 :     GEN L, S = gprec_w(s,LOWDEFAULTPREC), rS = real_i(S), iS = imag_i(S);
    1717             :     double l;
    1718             : 
    1719           7 :     l = rtodbl( gnorm(glog(S, 3)) );
    1720           7 :     l = log(l) / 2.;
    1721           7 :     lim = 2 + (long)ceil((prec2nbits_mul(prec, LOG2) - l) / (2*(1+log((double)la))));
    1722           7 :     if (lim < 2) lim = 2;
    1723             : 
    1724           7 :     l = (2*lim-1)*la / (2.*M_PI);
    1725           7 :     L = gsub(dbltor(l*l), gsqr(iS));
    1726           7 :     if (signe(L) < 0) L = gen_0;
    1727             : 
    1728           7 :     L = gsub(gsqrt(L, 3), rS);
    1729           7 :     if (signe(L) > 0) nn = (long)ceil(rtodbl(L)); else nn = 1;
    1730           7 :     if (DEBUGLEVEL>2) err_printf("lim, nn: [%ld, %ld]\n",lim,nn);
    1731             :   }
    1732             :   else
    1733             :   {
    1734          35 :     double ssig = rtodbl(sig);
    1735          35 :     double st = typ(s) == t_REAL? 0.0: rtodbl(imag_i(s));
    1736             :     double l;
    1737             :     {
    1738             :       double rlog, ilog; /* log (s - Euler) */
    1739          35 :       dcxlog(ssig - 0.57721566, st, &rlog,&ilog);
    1740          35 :       l = dnorm(rlog,ilog);
    1741             :     }
    1742          35 :     if (l < 0.000001) l = 0.000001;
    1743          35 :     l = log(l) / 2.;
    1744          35 :     lim = 2 + (long)ceil((prec2nbits_mul(prec, LOG2) - l) / (2*(1+log((double)la))));
    1745          35 :     if (lim < 2) lim = 2;
    1746             : 
    1747          35 :     l = (2*lim-1)*la / (2.*M_PI);
    1748          35 :     l = l*l - st*st;
    1749          35 :     if (l < 0.) l = 0.;
    1750          35 :     nn = (long)ceil( sqrt(l) - ssig );
    1751          35 :     if (nn < 1) nn = 1;
    1752          35 :     if (DEBUGLEVEL>2) err_printf("lim, nn: [%ld, %ld]\n",lim,nn);
    1753             :   }
    1754          42 :   incrprec(prec); unr = real_1(prec); /* one extra word of precision */
    1755          42 :   s2 = gmul2n(s, 1); sq = gsqr(s);
    1756          42 :   a = gdiv(unr, gaddgs(s, nn)); /* 1 / (s+n) */
    1757          42 :   av2 = avma; sum = gmul2n(a, -1);
    1758        1050 :   for (k = 0; k < nn - 1; k += 2)
    1759             :   {
    1760        1008 :     GEN tmp = gaddsg(k*(k + 1), gadd(gmulsg(2*k + 1, s), sq));
    1761        1008 :     sum = gadd(sum, gdiv(gaddsg(2*k + 1, s2), tmp));
    1762        1008 :     if ((k & 1023) == 0) sum = gerepileupto(av2, sum);
    1763             :   }
    1764          42 :   if (odd(nn)) sum = gadd(sum, gdiv(unr, gaddsg(nn - 1, s)));
    1765          42 :   z = gsub(glog(gaddgs(s, nn), prec), sum);
    1766          42 :   if (DEBUGLEVEL>2) timer_printf(&T,"sum from 0 to N - 1");
    1767             : 
    1768          42 :   in2 = gsqr(a);
    1769          42 :   mpbern(lim,prec);
    1770          42 :   av2 = avma; tes = divru(bernreal(2*lim, prec), 2*lim);
    1771        2317 :   for (k=2*lim-2; k>=2; k-=2)
    1772             :   {
    1773        2275 :     tes = gadd(gmul(in2,tes), divru(bernreal(k, prec), k));
    1774        2275 :     if ((k & 255) == 0) tes = gerepileupto(av2, tes);
    1775             :   }
    1776          42 :   if (DEBUGLEVEL>2) timer_printf(&T,"Bernoulli sum");
    1777          42 :   z = gsub(z, gmul(in2,tes));
    1778          42 :   if (funeq)
    1779             :   {
    1780          14 :     GEN pi = mppi(prec);
    1781          14 :     z = gadd(z, gmul(pi, gcotan(gmul(pi,s), prec)));
    1782             :   }
    1783          42 :   avma = av; return affc_fixlg(z, res);
    1784             : }
    1785             : 
    1786             : /* n > 0; return psi(1+x) + O(x^n), x = pol_x(v) */
    1787             : static GEN
    1788        6013 : serpsi1(long n, long v, long prec)
    1789             : {
    1790        6013 :   long i, l = n+3;
    1791        6013 :   GEN z, g, s = cgetg(l, t_SER);
    1792        6013 :   s[1] = evalsigne(1)|evalvalp(0)|evalvarn(v);
    1793        6013 :   g = mpeuler(prec); setsigne(g, -1);
    1794        6013 :   z = veczeta(gen_1, gen_2, n, prec); /* zeta(2..n) */
    1795        6013 :   gel(s,2) = g;
    1796       25417 :   for (i = 2; i < l-1; i++)
    1797             :   {
    1798       19404 :     GEN c = gel(z,i-1); /* zeta(i) */
    1799       19404 :     if (odd(i)) setsigne(c, -1);
    1800       19404 :     gel(s,i+1) = c;
    1801             :   }
    1802        6013 :   return s;
    1803             : }
    1804             : /* T an RgX, return T(X + z0) + O(X^L) */
    1805             : static GEN
    1806      788146 : tr(GEN T, GEN z0, long L)
    1807             : {
    1808      788146 :   GEN s = RgX_to_ser(RgX_translate(T, z0), L+3);
    1809      788146 :   setvarn(s, 0); return s;
    1810             : }
    1811             : /* z0 a complex number with Re(z0) > 1/2; return psi(z0+x) + O(x^L)
    1812             :  * using Luke's rational approximation for psi(x) */
    1813             : static GEN
    1814        3514 : serpsiz0(GEN z0, long L, long v, long prec)
    1815             : {
    1816             :   pari_sp av;
    1817             :   GEN A,A1,A2, B,B1,B2, Q;
    1818             :   long n;
    1819        3514 :   n = gprecision(z0); if (n) prec = n;
    1820        3514 :   z0 = gtofp(z0, prec + EXTRAPRECWORD);
    1821             :   /* Start from n = 3; in Luke's notation, A2 := A_{n-2}, A1 := A_{n-1},
    1822             :    * A := A_n. Same for B */
    1823        3514 :   av = avma;
    1824        3514 :   A2= gdivgs(mkpoln(2, gen_1, utoipos(6)), 2);
    1825        3514 :   B2 = scalarpol_shallow(utoipos(4), 0);
    1826        3514 :   A1= gdivgs(mkpoln(3, gen_1, utoipos(82), utoipos(96)), 6);
    1827        3514 :   B1 = mkpoln(2, utoipos(8), utoipos(28));
    1828        3514 :   A = gdivgs(mkpoln(4, gen_1, utoipos(387), utoipos(2906), utoipos(1920)), 12);
    1829        3514 :   B = mkpoln(3, utoipos(14), utoipos(204), utoipos(310));
    1830        3514 :   A2= tr(A2,z0, L);
    1831        3514 :   B2= tr(B2,z0, L);
    1832        3514 :   A1= tr(A1,z0, L);
    1833        3514 :   B1= tr(B1,z0, L);
    1834        3514 :   A = tr(A, z0, L);
    1835        3514 :   B = tr(B, z0, L); Q = gdiv(A, B);
    1836             :   /* work with z0+x as a variable */
    1837      254516 :   for (n = 4;; n++)
    1838             :   {
    1839      254516 :     GEN Q0 = Q, a, b, r, c3,c2,c1,c0 = muluu(2*n-3, n+1);
    1840      254516 :     GEN u = subiu(muluu(n, 7*n-9), 6);
    1841      254516 :     GEN t = addiu(muluu(n, 7*n-19), 4);
    1842             :     /* c1=(2*n-1)*(3*(n-1)*z+7*n^2-9*n-6);
    1843             :      * c2=(2*n-3)*(z-n-1)*(-3*(n-1)*z+7*n^2-19*n+4);
    1844             :      * c3=(2*n-1)*(n-3)*(z-n)*(z-(n+1))*(z+(n-4)); */
    1845      254516 :     c1 = deg1pol_shallow(muluu(3*(n-1),2*n-1), muliu(u,2*n-1), 0);
    1846      254516 :     c2 = ZX_mul(deg1pol_shallow(utoipos(2*n-3), negi(muluu(2*n-3,n+1)), 0),
    1847      254516 :                 deg1pol_shallow(utoineg(3*(n-1)), t, 0));
    1848      254516 :     r = mkvec3(utoipos(n), utoipos(n+1), stoi(4-n));
    1849      254516 :     c3 = ZX_Z_mul(roots_to_pol(r,0), muluu(2*n-1,n-3));
    1850      254516 :     c1 = tr(c1, z0, L+3);
    1851      254516 :     c2 = tr(c2, z0, L+3);
    1852      254516 :     c3 = tr(c3, z0, L+3);
    1853             : 
    1854             :     /* A_{n+1}, B_{n+1} */
    1855      254516 :     a = gdiv(gadd(gadd(gmul(c1,A),gmul(c2,A1)),gmul(c3,A2)), c0);
    1856      254516 :     b = gdiv(gadd(gadd(gmul(c1,B),gmul(c2,B1)),gmul(c3,B2)), c0);
    1857      254516 :     Q = gdiv(a,b);
    1858      254516 :     if (gexpo(gsub(Q,Q0)) < -prec2nbits(prec)) break;
    1859      251002 :     A2 = A1; A1 = A; A = a;
    1860      251002 :     B2 = B1; B1 = B; B = b;
    1861      251002 :     if (gc_needed(av,1))
    1862             :     {
    1863           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"serpsiz0, n = %ld", n);
    1864           0 :       gerepileall(av, 7, &A,&A1,&A2, &B,&B1,&B2, &Q);
    1865             :     }
    1866      251002 :   }
    1867        3514 :   Q = gmul(Q, gmul2n(gsubsg(1, ginv(tr(pol_x(v),z0, L))), 1));
    1868        3514 :   setvarn(Q, v);
    1869        3514 :   return gadd(negr(mpeuler(prec)), Q);
    1870             : }
    1871             : /* sum (-1)^k*H(m,k)x^k + O(x^L); L > 0;
    1872             :  * H(m,k) = (-1)^{k * \delta_{m > 0}} sum_{1<=i<m} 1/i^(k+1) */
    1873             : static GEN
    1874         511 : Hseries(long m, long L, long v, long prec)
    1875             : {
    1876         511 :   long i, k, bit, l = L+3, M = m < 0? 1-m: m;
    1877         511 :   pari_sp av = avma;
    1878         511 :   GEN H = cgetg(l, t_SER);
    1879         511 :   H[1] = evalsigne(1)|evalvarn(v)|evalvalp(0);
    1880         511 :   prec++;
    1881         511 :   bit = -prec2nbits(prec);
    1882         511 :   for(k = 2; k < l; k++) gel(H,k) = gen_1; /* i=1 */
    1883         532 :   for (i = 2; i < M; i++)
    1884             :   {
    1885          21 :     GEN ik = invr(utor(i, prec));
    1886         154 :     for (k = 2; k < l; k++)
    1887             :     {
    1888         133 :       if (k > 2) { ik = divru(ik, i); if (expo(ik) < bit) break; }
    1889         133 :       gel(H,k) = gadd(gel(H,k), ik);
    1890             :     }
    1891          21 :     if (gc_needed(av,3))
    1892             :     {
    1893           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"Hseries, i = %ld/%ld", i,M);
    1894           0 :       H = gerepilecopy(av, H);
    1895             :     }
    1896             :   }
    1897         511 :   if (m > 0)
    1898         497 :     for (k = 3; k < l; k+=2) togglesign_safe(&gel(H,k));
    1899         511 :   return H;
    1900             : }
    1901             : 
    1902             : static GEN
    1903        9541 : serpsi(GEN y, long prec)
    1904             : {
    1905        9541 :   GEN Q = NULL, z0, Y = y, Y2;
    1906        9541 :   long L = lg(y)-2, v  = varn(y), vy = valp(y);
    1907             : 
    1908        9541 :   if (!L) pari_err_DOMAIN("psi", "argument", "=", gen_0,y);
    1909        9534 :   if (vy < 0) pari_err_DOMAIN("psi", "series valuation", "<", gen_0,y);
    1910        9534 :   if (vy)
    1911          14 :     z0 = gen_0;
    1912             :   else
    1913             :   {
    1914        9520 :     z0 = simplify_shallow(gel(y,2));
    1915        9520 :     (void)isint(z0, &z0);
    1916             :   }
    1917        9534 :   if (typ(z0) == t_INT && !is_bigint(z0))
    1918             :   {
    1919        6020 :     long m = itos(z0);
    1920        6020 :     if (abscmpiu(muluu(prec2nbits(prec),L), labs(m)) > 0)
    1921             :     { /* psi(m+x) = psi(1+x) + sum_{1 <= i < m} 1/(i+x) for m > 0
    1922             :                     psi(1+x) - sum_{0 <= i < -m} 1/(i+x) for m <= 0 */
    1923        6020 :       GEN H = NULL;
    1924        6020 :       if (m <= 0) L--; /* lose series accuracy due to 1/x term */
    1925        6020 :       if (L)
    1926             :       {
    1927        6013 :         Q = serpsi1(L, v, prec);
    1928        6013 :         if (m && m != 1) { H = Hseries(m, L, v, prec); Q = gadd(Q, H); }
    1929        6013 :         if (m <= 0) Q = gsub(Q, ginv(pol_x(v)));
    1930             :       }
    1931             :       else
    1932             :       {
    1933           7 :         Q = scalarser(gen_m1, v, 1);
    1934           7 :         setvalp(Q,-1);
    1935             :       }
    1936             :     }
    1937             :   }
    1938        9534 :   if (!Q)
    1939             :   { /* use psi(1-y)=psi(y)+Pi*cotan(Pi*y) ? */
    1940        3514 :     if (gcmp(real_i(z0),ghalf) < 0) { z0 = gsubsg(1,z0); Y = gsubsg(1,y); }
    1941        3514 :     Q = serpsiz0(z0, L, v, prec);
    1942             :   }
    1943        9534 :   Y2 = serchop0(Y); if (signe(Y2)) Q = gsubst(Q, v, Y2);
    1944             :   /* psi(z0 + Y2) = psi(Y) */
    1945        9534 :   if (Y != y)
    1946             :   { /* psi(y) = psi(Y) + Pi cotan(Pi Y) */
    1947          49 :     GEN pi = mppi(prec);
    1948          49 :     if (typ(z0) == t_INT) Y = Y2; /* in this case cotan(Pi*Y2) = cotan(Pi*Y) */
    1949          49 :     Q = gadd(Q, gmul(pi, gcotan(gmul(pi,Y), prec)));
    1950             :   }
    1951        9534 :   return Q;
    1952             : }
    1953             : 
    1954             : GEN
    1955        9597 : gpsi(GEN x, long prec)
    1956             : {
    1957             :   pari_sp av;
    1958             :   GEN y;
    1959        9597 :   switch(typ(x))
    1960             :   {
    1961          42 :     case t_REAL: case t_COMPLEX: return cxpsi(x,prec);
    1962             :     default:
    1963        9555 :       av = avma; if (!(y = toser_i(x))) break;
    1964        9541 :       return gerepileupto(av, serpsi(y,prec));
    1965             :   }
    1966          14 :   return trans_eval("psi",gpsi,x,prec);
    1967             : }

Generated by: LCOV version 1.11