Code coverage tests

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

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

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

LCOV - code coverage report
Current view: top level - basemath - trans2.c (source / functions) Hit Total Coverage
Test: PARI/GP v2.12.1 lcov report (development 25819-e703fe1174) Lines: 1140 1192 95.6 %
Date: 2020-09-18 06:10:04 Functions: 87 89 97.8 %
Legend: Lines: hit not hit

          Line data    Source code
       1             : /* Copyright (C) 2000  The PARI group.
       2             : 
       3             : This file is part of the PARI/GP package.
       4             : 
       5             : PARI/GP is free software; you can redistribute it and/or modify it under the
       6             : terms of the GNU General Public License as published by the Free Software
       7             : Foundation. It is distributed in the hope that it will be useful, but WITHOUT
       8             : ANY WARRANTY WHATSOEVER.
       9             : 
      10             : Check the License for details. You should have received a copy of it, along
      11             : with the package; see the file 'COPYING'. If not, write to the Free Software
      12             : Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. */
      13             : 
      14             : /********************************************************************/
      15             : /**                                                                **/
      16             : /**                   TRANSCENDENTAL FUNCTIONS                     **/
      17             : /**                          (part 2)                              **/
      18             : /**                                                                **/
      19             : /********************************************************************/
      20             : #include "pari.h"
      21             : #include "paripriv.h"
      22             : 
      23             : GEN
      24      244487 : trans_fix_arg(long *prec, GEN *s0, GEN *sig, GEN *tau, pari_sp *av, GEN *res)
      25             : {
      26      244487 :   GEN p1, s = *s0 = cxtoreal(*s0);
      27             :   long l;
      28      244487 :   l = precision(s); if (!l) l = *prec;
      29      244487 :   if (l < LOWDEFAULTPREC) l = LOWDEFAULTPREC;
      30      244487 :   *res = cgetc(l); *av = avma;
      31      244492 :   if (typ(s) == t_COMPLEX)
      32             :   { /* s = sig + i t */
      33      210553 :     s = cxtofp(s, l+EXTRAPRECWORD);
      34      210553 :     *sig = gel(s,1);
      35      210553 :     *tau = gel(s,2);
      36             :   }
      37             :   else /* real number */
      38             :   {
      39       33939 :     *sig = s = gtofp(s, l+EXTRAPRECWORD);
      40       33933 :     *tau = gen_0;
      41       33933 :     p1 = trunc2nr(s, 0);
      42       33931 :     if (!signe(subri(s,p1))) *s0 = p1;
      43             :   }
      44      244478 :   *prec = l; return s;
      45             : }
      46             : 
      47             : /********************************************************************/
      48             : /**                                                                **/
      49             : /**                          ARCTANGENT                            **/
      50             : /**                                                                **/
      51             : /********************************************************************/
      52             : /* atan(b/a), real a and b, suitable for gerepileupto */
      53             : static GEN
      54         196 : atan2_agm(GEN a, GEN b, long prec)
      55         196 : { return gel(logagmcx(mkcomplex(a, b), prec), 2); }
      56             : static GEN
      57      722059 : mpatan(GEN x)
      58             : {
      59      722059 :   long l, l1, l2, n, m, i, lp, e, s, sx = signe(x);
      60             :   pari_sp av0, av;
      61             :   double alpha, beta, delta;
      62             :   GEN y, p1, p2, p3, p4, p5, unr;
      63             :   int inv;
      64             : 
      65      722059 :   if (!sx) return real_0_bit(expo(x));
      66      722024 :   l = lp = realprec(x);
      67      722024 :   if (absrnz_equal1(x)) { /* |x| = 1 */
      68        3730 :     y = Pi2n(-2, l+EXTRAPRECWORD); if (sx < 0) setsigne(y,-1);
      69        3730 :     return y;
      70             :   }
      71      718294 :   if (l > AGM_ATAN_LIMIT)
      72         175 :   { av = avma; return gerepileuptoleaf(av, atan2_agm(gen_1, x, l)); }
      73             : 
      74      718119 :   e = expo(x); inv = (e >= 0); /* = (|x| > 1 ) */
      75      718119 :   if (e > 0) lp += nbits2extraprec(e);
      76             : 
      77      718119 :   y = cgetr(lp); av0 = avma;
      78      718119 :   p1 = rtor(x, l+EXTRAPRECWORD); setabssign(p1); /* p1 = |x| */
      79      718119 :   if (inv) p1 = invr(p1);
      80      718119 :   e = expo(p1);
      81      718119 :   if (e < -100)
      82       10308 :     alpha = 1.65149612947 - e; /* log_2(Pi) - e */
      83             :   else
      84      707811 :     alpha = log2(M_PI / atan(rtodbl(p1)));
      85      718119 :   beta = (double)(prec2nbits(l)>>1);
      86      718119 :   delta = 1 + beta - alpha/2;
      87      718119 :   if (delta <= 0) { n = 1; m = 0; }
      88             :   else
      89             :   {
      90      710583 :     double fi = alpha-2;
      91      710583 :     if (delta >= fi*fi)
      92             :     {
      93      651328 :       double t = 1 + sqrt(delta);
      94      651328 :       n = (long)t;
      95      651328 :       m = (long)(t - fi);
      96             :     }
      97             :     else
      98             :     {
      99       59255 :       n = (long)(1+beta/fi);
     100       59255 :       m = 0;
     101             :     }
     102             :   }
     103      718119 :   l2 = l + nbits2extraprec(m);
     104      718119 :   p2 = rtor(p1, l2); av = avma;
     105     6981636 :   for (i=1; i<=m; i++)
     106             :   {
     107     6263517 :     p5 = addsr(1, sqrr(p2)); setprec(p5,l2);
     108     6263517 :     p5 = addsr(1, sqrtr_abs(p5)); setprec(p5,l2);
     109     6263517 :     affrr(divrr(p2,p5), p2); set_avma(av);
     110             :   }
     111      718119 :   p3 = sqrr(p2); l1 = minss(LOWDEFAULTPREC+EXTRAPRECWORD, l2); /* l1 increases to l2 */;
     112      718119 :   unr = real_1(l2); setprec(unr,l1);
     113      718119 :   p4 = cgetr(l2); setprec(p4,l1);
     114      718119 :   affrr(divru(unr,2*n+1), p4);
     115      718119 :   s = 0; e = expo(p3); av = avma;
     116     7909045 :   for (i = n; i > 1; i--) /* n >= 1. i = 1 done outside for efficiency */
     117             :   {
     118     7190926 :     setprec(p3,l1); p5 = mulrr(p4,p3);
     119     7190926 :     l1 += dvmdsBIL(s - e, &s); if (l1 > l2) l1 = l2;
     120     7190926 :     setprec(unr,l1); p5 = subrr(divru(unr,2*i-1), p5);
     121     7190926 :     setprec(p4,l1); affrr(p5,p4); set_avma(av);
     122             :   }
     123      718119 :   setprec(p3, l2); p5 = mulrr(p4,p3); /* i = 1 */
     124      718119 :   setprec(unr,l2); p4 = subrr(unr, p5);
     125             : 
     126      718119 :   p4 = mulrr(p2,p4); shiftr_inplace(p4, m);
     127      718119 :   if (inv) p4 = subrr(Pi2n(-1, lp), p4);
     128      718119 :   if (sx < 0) togglesign(p4);
     129      718119 :   affrr_fixlg(p4,y); set_avma(av0); return y;
     130             : }
     131             : 
     132             : GEN
     133       21355 : gatan(GEN x, long prec)
     134             : {
     135             :   pari_sp av;
     136             :   GEN a, y;
     137             : 
     138       21355 :   switch(typ(x))
     139             :   {
     140       12927 :     case t_REAL: return mpatan(x);
     141        7672 :     case t_COMPLEX: /* atan(x) = -i atanh(ix) */
     142        7672 :       if (ismpzero(gel(x,2))) return gatan(gel(x,1), prec);
     143        7441 :       av = avma; return gerepilecopy(av, mulcxmI(gatanh(mulcxI(x),prec)));
     144         756 :     default:
     145         756 :       av = avma; if (!(y = toser_i(x))) break;
     146          28 :       if (valp(y) < 0) pari_err_DOMAIN("atan","valuation", "<", gen_0, x);
     147          21 :       if (lg(y)==2) return gerepilecopy(av, y);
     148             :       /* lg(y) > 2 */
     149          14 :       a = integser(gdiv(derivser(y), gaddsg(1,gsqr(y))));
     150          14 :       if (!valp(y)) a = gadd(a, gatan(gel(y,2),prec));
     151          14 :       return gerepileupto(av, a);
     152             :   }
     153         728 :   return trans_eval("atan",gatan,x,prec);
     154             : }
     155             : /********************************************************************/
     156             : /**                                                                **/
     157             : /**                             ARCSINE                            **/
     158             : /**                                                                **/
     159             : /********************************************************************/
     160             : /* |x| < 1, x != 0 */
     161             : static GEN
     162          98 : mpasin(GEN x) {
     163          98 :   pari_sp av = avma;
     164          98 :   GEN z, a = sqrtr(subsr(1, sqrr(x)));
     165          98 :   if (realprec(x) > AGM_ATAN_LIMIT)
     166           7 :     z = atan2_agm(a, x, realprec(x));
     167             :   else
     168          91 :     z = mpatan(divrr(x, a));
     169          98 :   return gerepileuptoleaf(av, z);
     170             : }
     171             : 
     172             : static GEN mpacosh(GEN x);
     173             : GEN
     174        8183 : gasin(GEN x, long prec)
     175             : {
     176             :   long sx;
     177             :   pari_sp av;
     178             :   GEN a, y, p1;
     179             : 
     180        8183 :   switch(typ(x))
     181             :   {
     182         483 :     case t_REAL: sx = signe(x);
     183         483 :       if (!sx) return real_0_bit(expo(x));
     184         476 :       if (absrnz_equal1(x)) { /* |x| = 1 */
     185          28 :         if (sx > 0) return Pi2n(-1, realprec(x)); /* 1 */
     186          14 :         y = Pi2n(-1, realprec(x)); setsigne(y, -1); return y; /* -1 */
     187             :       }
     188         448 :       if (expo(x) < 0) return mpasin(x);
     189         350 :       y = cgetg(3,t_COMPLEX);
     190         350 :       gel(y,1) = Pi2n(-1, realprec(x));
     191         350 :       gel(y,2) = mpacosh(x);
     192         350 :       if (sx < 0) togglesign(gel(y,1)); else togglesign(gel(y,2));
     193         350 :       return y;
     194             : 
     195        7637 :     case t_COMPLEX: /* asin(z) = -i asinh(iz) */
     196        7637 :       if (ismpzero(gel(x,2))) return gasin(gel(x,1), prec);
     197        7406 :       av = avma;
     198        7406 :       return gerepilecopy(av, mulcxmI(gasinh(mulcxI(x), prec)));
     199          63 :     default:
     200          63 :       av = avma; if (!(y = toser_i(x))) break;
     201          42 :       if (gequal0(y)) return gerepilecopy(av, y);
     202             :       /* lg(y) > 2*/
     203          35 :       if (valp(y) < 0) pari_err_DOMAIN("asin","valuation", "<", gen_0, x);
     204          28 :       p1 = gsubsg(1,gsqr(y));
     205          28 :       if (gequal0(p1))
     206             :       {
     207          21 :         GEN t = Pi2n(-1,prec);
     208          21 :         if (gsigne(gel(y,2)) < 0) setsigne(t, -1);
     209          21 :         return gerepileupto(av, scalarser(t, varn(y), valp(p1)>>1));
     210             :       }
     211           7 :       p1 = gdiv(derivser(y), gsqrt(p1,prec));
     212           7 :       a = integser(p1);
     213           7 :       if (!valp(y)) a = gadd(a, gasin(gel(y,2),prec));
     214           7 :       return gerepileupto(av, a);
     215             :   }
     216          21 :   return trans_eval("asin",gasin,x,prec);
     217             : }
     218             : /********************************************************************/
     219             : /**                                                                **/
     220             : /**                             ARCCOSINE                          **/
     221             : /**                                                                **/
     222             : /********************************************************************/
     223             : static GEN
     224          14 : acos0(long e) { return Pi2n(-1, nbits2prec(e<0? -e: 1)); }
     225             : 
     226             : /* |x| < 1, x != 0 */
     227             : static GEN
     228         105 : mpacos(GEN x)
     229             : {
     230         105 :   pari_sp av = avma;
     231         105 :   GEN z, a = sqrtr(subsr(1, sqrr(x)));
     232         105 :   if (realprec(x) > AGM_ATAN_LIMIT)
     233          14 :     z = atan2_agm(x, a, realprec(x));
     234             :   else
     235             :   {
     236          91 :     z = mpatan(divrr(a, x));
     237          91 :     if (signe(x) < 0) z = addrr(mppi(realprec(z)), z);
     238             :   }
     239         105 :   return gerepileuptoleaf(av, z);
     240             : }
     241             : 
     242             : GEN
     243        7938 : gacos(GEN x, long prec)
     244             : {
     245             :   long sx;
     246             :   pari_sp av;
     247             :   GEN a, y, p1;
     248             : 
     249        7938 :   switch(typ(x))
     250             :   {
     251         252 :     case t_REAL: sx = signe(x);
     252         252 :       if (!sx) return acos0(expo(x));
     253         245 :       if (absrnz_equal1(x)) /* |x| = 1 */
     254          14 :         return sx > 0? real_0_bit( -(bit_prec(x)>>1) ) : mppi(realprec(x));
     255         231 :       if (expo(x) < 0) return mpacos(x);
     256             : 
     257         175 :       y = cgetg(3,t_COMPLEX); p1 = mpacosh(x);
     258         175 :       if (sx < 0) { gel(y,1) = mppi(realprec(x)); togglesign(p1); }
     259          91 :       else gel(y,1) = gen_0;
     260         175 :       gel(y,2) = p1; return y;
     261             : 
     262        7637 :     case t_COMPLEX:
     263        7637 :       if (ismpzero(gel(x,2))) return gacos(gel(x,1), prec);
     264        7406 :       av = avma;
     265        7406 :       p1 = gadd(x, mulcxI(gsqrt(gsubsg(1,gsqr(x)), prec)));
     266        7406 :       y = glog(p1,prec); /* log(x + I*sqrt(1-x^2)) */
     267        7406 :       return gerepilecopy(av, mulcxmI(y));
     268          49 :     default:
     269          49 :       av = avma; if (!(y = toser_i(x))) break;
     270          35 :       if (valp(y) < 0) pari_err_DOMAIN("acos","valuation", "<", gen_0, x);
     271          28 :       if (lg(y) > 2)
     272             :       {
     273          21 :         p1 = gsubsg(1,gsqr(y));
     274          21 :         if (gequal0(p1)) { set_avma(av); return zeroser(varn(y), valp(p1)>>1); }
     275           7 :         p1 = integser(gdiv(gneg(derivser(y)), gsqrt(p1,prec)));
     276             :         /*y(t) = 1+O(t)*/
     277           7 :         if (gequal1(gel(y,2)) && !valp(y)) return gerepileupto(av, p1);
     278             :       }
     279           7 :       else p1 = y;
     280          14 :       a = (lg(y)==2 || valp(y))? Pi2n(-1, prec): gacos(gel(y,2),prec);
     281          14 :       return gerepileupto(av, gadd(a,p1));
     282             :   }
     283          14 :   return trans_eval("acos",gacos,x,prec);
     284             : }
     285             : /********************************************************************/
     286             : /**                                                                **/
     287             : /**                            ARGUMENT                            **/
     288             : /**                                                                **/
     289             : /********************************************************************/
     290             : 
     291             : /* we know that x and y are not both 0 */
     292             : static GEN
     293      708970 : mparg(GEN x, GEN y)
     294             : {
     295      708970 :   long prec, sx = signe(x), sy = signe(y);
     296             :   GEN z;
     297             : 
     298      708970 :   if (!sy)
     299             :   {
     300           0 :     if (sx > 0) return real_0_bit(expo(y) - expo(x));
     301           0 :     return mppi(realprec(x));
     302             :   }
     303      708970 :   prec = realprec(y); if (prec < realprec(x)) prec = realprec(x);
     304      708970 :   if (!sx)
     305             :   {
     306          20 :     z = Pi2n(-1, prec); if (sy < 0) setsigne(z,-1);
     307          20 :     return z;
     308             :   }
     309             : 
     310      708950 :   if (expo(x)-expo(y) > -2)
     311             :   {
     312      621662 :     z = mpatan(divrr(y,x)); if (sx > 0) return z;
     313      137804 :     return addrr_sign(z, signe(z), mppi(prec), sy);
     314             :   }
     315       87288 :   z = mpatan(divrr(x,y));
     316       87288 :   return addrr_sign(z, -signe(z), Pi2n(-1, prec), sy);
     317             : }
     318             : 
     319             : static GEN
     320     1417940 : rfix(GEN x,long prec)
     321             : {
     322     1417940 :   switch(typ(x))
     323             :   {
     324        3139 :     case t_INT: return itor(x, prec);
     325       39032 :     case t_FRAC: return fractor(x, prec);
     326     1375769 :     case t_REAL: break;
     327           0 :     default: pari_err_TYPE("rfix (conversion to t_REAL)",x);
     328             :   }
     329     1375769 :   return x;
     330             : }
     331             : 
     332             : static GEN
     333      708970 : cxarg(GEN x, GEN y, long prec)
     334             : {
     335      708970 :   pari_sp av = avma;
     336      708970 :   x = rfix(x,prec);
     337      708970 :   y = rfix(y,prec); return gerepileuptoleaf(av, mparg(x,y));
     338             : }
     339             : 
     340             : GEN
     341      717370 : garg(GEN x, long prec)
     342             : {
     343             :   long l;
     344      717370 :   if (gequal0(x)) pari_err_DOMAIN("arg", "argument", "=", gen_0, x);
     345      717370 :   switch(typ(x))
     346             :   {
     347        8400 :     case t_REAL: prec = realprec(x); /* fall through */
     348        8400 :     case t_INT: case t_FRAC: return (gsigne(x)>0)? real_0(prec): mppi(prec);
     349      708970 :     case t_COMPLEX:
     350      708970 :       l = precision(x); if (l) prec = l;
     351      708970 :       return cxarg(gel(x,1),gel(x,2),prec);
     352             :   }
     353           0 :   return trans_eval("arg",garg,x,prec);
     354             : }
     355             : 
     356             : /********************************************************************/
     357             : /**                                                                **/
     358             : /**                      HYPERBOLIC COSINE                         **/
     359             : /**                                                                **/
     360             : /********************************************************************/
     361             : /* 1 + x */
     362             : static GEN
     363           7 : mpcosh0(long e) { return e >= 0? real_0_bit(e): real_1_bit(-e); }
     364             : static GEN
     365        3696 : mpcosh(GEN x)
     366             : {
     367             :   pari_sp av;
     368             :   GEN z;
     369             : 
     370        3696 :   if (!signe(x)) return mpcosh0(expo(x));
     371        3689 :   av = avma;
     372        3689 :   z = mpexp(x); z = addrr(z, invr(z)); shiftr_inplace(z, -1);
     373        3689 :   return gerepileuptoleaf(av, z);
     374             : }
     375             : 
     376             : GEN
     377        3780 : gcosh(GEN x, long prec)
     378             : {
     379             :   pari_sp av;
     380             :   GEN y, p1;
     381             : 
     382        3780 :   switch(typ(x))
     383             :   {
     384        3696 :     case t_REAL: return mpcosh(x);
     385          21 :     case t_COMPLEX:
     386          21 :       if (isintzero(gel(x,1))) return gcos(gel(x,2),prec);
     387             :       /* fall through */
     388             :     case t_PADIC:
     389          21 :       av = avma; p1 = gexp(x,prec); p1 = gadd(p1, ginv(p1));
     390          21 :       return gerepileupto(av, gmul2n(p1,-1));
     391          49 :     default:
     392          49 :       av = avma; if (!(y = toser_i(x))) break;
     393          28 :       if (gequal0(y) && valp(y) == 0) return gerepilecopy(av, y);
     394          28 :       p1 = gexp(y,prec); p1 = gadd(p1, ginv(p1));
     395          28 :       return gerepileupto(av, gmul2n(p1,-1));
     396             :   }
     397          21 :   return trans_eval("cosh",gcosh,x,prec);
     398             : }
     399             : /********************************************************************/
     400             : /**                                                                **/
     401             : /**                       HYPERBOLIC SINE                          **/
     402             : /**                                                                **/
     403             : /********************************************************************/
     404             : static GEN
     405           0 : mpsinh0(long e) { return real_0_bit(e); }
     406             : static GEN
     407       15351 : mpsinh(GEN x)
     408             : {
     409             :   pari_sp av;
     410       15351 :   long ex = expo(x), lx;
     411             :   GEN z, res;
     412             : 
     413       15351 :   if (!signe(x)) return mpsinh0(ex);
     414       15351 :   lx = realprec(x); res = cgetr(lx); av = avma;
     415       15351 :   if (ex < 1 - BITS_IN_LONG)
     416             :   { /* y = e^x-1; e^x - e^(-x) = y(1 + 1/(y+1)) */
     417           7 :     GEN y = mpexpm1(x);
     418           7 :     z = addrs(y,1); if (realprec(z) > lx+1) z = rtor(z,lx+1); /* e^x */
     419           7 :     z = mulrr(y, addsr(1,invr(z)));
     420             :   }
     421             :   else
     422             :   {
     423       15344 :     z = mpexp(x);
     424       15344 :     z = subrr(z, invr(z));
     425             :   }
     426       15351 :   shiftr_inplace(z, -1);
     427       15351 :   affrr(z, res); set_avma(av); return res;
     428             : }
     429             : 
     430             : void
     431       17164 : mpsinhcosh(GEN x, GEN *s, GEN *c)
     432             : {
     433             :   pari_sp av;
     434       17164 :   long lx, ex = expo(x);
     435             :   GEN z, zi, S, C;
     436       17164 :   if (!signe(x))
     437             :   {
     438           0 :     *c = mpcosh0(ex);
     439           0 :     *s = mpsinh0(ex); return;
     440             :   }
     441       17164 :   lx = realprec(x);
     442       17164 :   *c = cgetr(lx);
     443       17164 :   *s = cgetr(lx); av = avma;
     444       17164 :   if (ex < 1 - BITS_IN_LONG)
     445             :   { /* y = e^x-1; e^x - e^(-x) = y(1 + 1/(y+1)) */
     446         403 :     GEN y = mpexpm1(x);
     447         403 :     z = addrs(y,1); if (realprec(z) > lx+1) z = rtor(z,lx+1); /* e^x */
     448         403 :     zi = invr(z); /* z = exp(x), zi = exp(-x) */
     449         403 :     S = mulrr(y, addsr(1,zi));
     450             :   }
     451             :   else
     452             :   {
     453       16761 :     z = mpexp(x);
     454       16761 :     zi = invr(z);
     455       16761 :     S = subrr(z, zi);
     456             :   }
     457       17164 :   C = addrr(z, zi);
     458       17164 :   shiftr_inplace(S, -1); affrr(S, *s);
     459       17164 :   shiftr_inplace(C, -1); affrr(C, *c); set_avma(av);
     460             : }
     461             : 
     462             : GEN
     463       17773 : gsinh(GEN x, long prec)
     464             : {
     465             :   pari_sp av;
     466             :   GEN y, p1;
     467             : 
     468       17773 :   switch(typ(x))
     469             :   {
     470       15351 :     case t_REAL: return mpsinh(x);
     471          21 :     case t_COMPLEX:
     472          21 :       if (isintzero(gel(x,1))) retmkcomplex(gen_0, gsin(gel(x,2),prec));
     473             :       /* fall through */
     474             :     case t_PADIC:
     475          14 :       av = avma; p1 = gexp(x,prec); p1 = gsub(p1, ginv(p1));
     476          14 :       return gerepileupto(av, gmul2n(p1,-1));
     477        2394 :     default:
     478        2394 :       av = avma; if (!(y = toser_i(x))) break;
     479        2366 :       if (gequal0(y) && valp(y) == 0) return gerepilecopy(av, y);
     480        2366 :       p1 = gexp(y, prec); p1 = gsub(p1, ginv(p1));
     481        2366 :       return gerepileupto(av, gmul2n(p1,-1));
     482             :   }
     483          28 :   return trans_eval("sinh",gsinh,x,prec);
     484             : }
     485             : /********************************************************************/
     486             : /**                                                                **/
     487             : /**                      HYPERBOLIC TANGENT                        **/
     488             : /**                                                                **/
     489             : /********************************************************************/
     490             : 
     491             : static GEN
     492       77056 : mptanh(GEN x)
     493             : {
     494       77056 :   long lx, s = signe(x);
     495             :   GEN y;
     496             : 
     497       77056 :   if (!s) return real_0_bit(expo(x));
     498       77056 :   lx = realprec(x);
     499       77056 :   if (abscmprr(x, utor(prec2nbits(lx), LOWDEFAULTPREC)) >= 0) {
     500       24840 :     y = real_1(lx);
     501             :   } else {
     502       52216 :     pari_sp av = avma;
     503       52216 :     long ex = expo(x);
     504             :     GEN t;
     505       52216 :     if (ex < 1 - BITS_IN_LONG) x = rtor(x, lx + nbits2extraprec(-ex)-1);
     506       52216 :     t = exp1r_abs(gmul2n(x,1)); /* exp(|2x|) - 1 */
     507       52216 :     y = gerepileuptoleaf(av, divrr(t, addsr(2,t)));
     508             :   }
     509       77056 :   if (s < 0) togglesign(y); /* tanh is odd */
     510       77056 :   return y;
     511             : }
     512             : 
     513             : GEN
     514       77161 : gtanh(GEN x, long prec)
     515             : {
     516             :   pari_sp av;
     517             :   GEN y, t;
     518             : 
     519       77161 :   switch(typ(x))
     520             :   {
     521       77056 :     case t_REAL: return mptanh(x);
     522          35 :     case t_COMPLEX:
     523          35 :       if (isintzero(gel(x,1))) retmkcomplex(gen_0, gtan(gel(x,2),prec));
     524             :       /* fall through */
     525             :     case t_PADIC:
     526          28 :       av = avma;
     527          28 :       t = gexp(gmul2n(x,1),prec);
     528          28 :       t = gdivsg(-2, gaddgs(t,1));
     529          28 :       return gerepileupto(av, gaddsg(1,t));
     530          63 :     default:
     531          63 :       av = avma; if (!(y = toser_i(x))) break;
     532          28 :       if (gequal0(y)) return gerepilecopy(av, y);
     533          14 :       t = gexp(gmul2n(y, 1),prec);
     534          14 :       t = gdivsg(-2, gaddgs(t,1));
     535          14 :       return gerepileupto(av, gaddsg(1,t));
     536             :   }
     537          35 :   return trans_eval("tanh",gtanh,x,prec);
     538             : }
     539             : 
     540             : static GEN
     541           7 : mpcotanh(GEN x)
     542             : {
     543           7 :   long lx, s = signe(x);
     544             :   GEN y;
     545             : 
     546           7 :   if (!s) pari_err_DOMAIN("cotan", "argument", "=", gen_0, x);
     547             : 
     548           7 :   lx = realprec(x);
     549           7 :   if (abscmprr(x, utor(prec2nbits(lx), LOWDEFAULTPREC)) >= 0) {
     550           0 :     y = real_1(lx);
     551             :   } else {
     552           7 :     pari_sp av = avma;
     553           7 :     long ex = expo(x);
     554             :     GEN t;
     555           7 :     if (ex < 1 - BITS_IN_LONG) x = rtor(x, lx + nbits2extraprec(-ex)-1);
     556           7 :     t = exp1r_abs(gmul2n(x,1)); /* exp(|2x|) - 1 */
     557           7 :     y = gerepileuptoleaf(av, divrr(addsr(2,t), t));
     558             :   }
     559           7 :   if (s < 0) togglesign(y); /* cotanh is odd */
     560           7 :   return y;
     561             : }
     562             : 
     563             : GEN
     564          63 : gcotanh(GEN x, long prec)
     565             : {
     566             :   pari_sp av;
     567             :   GEN y, t;
     568             : 
     569          63 :   switch(typ(x))
     570             :   {
     571           7 :     case t_REAL: return mpcotanh(x);
     572          14 :     case t_COMPLEX:
     573          14 :       if (isintzero(gel(x,1))) retmkcomplex(gen_0, gcotan(gel(x,2),prec));
     574             :       /* fall through */
     575             :     case t_PADIC:
     576          14 :       av = avma;
     577          14 :       t = gexpm1(gmul2n(x,1),prec);
     578          14 :       return gerepileupto(av, gaddsg(1, gdivsg(2,t)));
     579          35 :     default:
     580          35 :       av = avma; if (!(y = toser_i(x))) break;
     581          28 :       if (gequal0(y)) return gerepilecopy(av, y);
     582          14 :       t = gexpm1(gmul2n(y,1),prec);
     583          14 :       return gerepileupto(av, gaddsg(1, gdivsg(2,t)));
     584             :   }
     585           7 :   return trans_eval("cotanh",gcotanh,x,prec);
     586             : }
     587             : 
     588             : /********************************************************************/
     589             : /**                                                                **/
     590             : /**                     AREA HYPERBOLIC SINE                       **/
     591             : /**                                                                **/
     592             : /********************************************************************/
     593             : 
     594             : /* x != 0 */
     595             : static GEN
     596         483 : mpasinh(GEN x)
     597             : {
     598         483 :   long lx = realprec(x), ex = expo(x);
     599         483 :   GEN z, res = cgetr(lx);
     600         483 :   pari_sp av = avma;
     601         483 :   if (ex < 1 - BITS_IN_LONG) x = rtor(x, lx + nbits2extraprec(-ex)-1);
     602         483 :   z = logr_abs( addrr_sign(x,1, sqrtr_abs( addrs(sqrr(x), 1) ), 1) );
     603         483 :   if (signe(x) < 0) togglesign(z);
     604         483 :   affrr(z, res); return gc_const(av, res);
     605             : }
     606             : 
     607             : GEN
     608       15729 : gasinh(GEN x, long prec)
     609             : {
     610             :   pari_sp av;
     611             :   GEN a, y, p1;
     612             : 
     613       15729 :   switch(typ(x))
     614             :   {
     615         490 :     case t_REAL:
     616         490 :       if (!signe(x)) return rcopy(x);
     617         483 :       return mpasinh(x);
     618             : 
     619       15043 :     case t_COMPLEX: {
     620             :       GEN a, b, d;
     621       15043 :       if (ismpzero(gel(x,2))) return gasinh(gel(x,1), prec);
     622       14588 :       av = avma;
     623       14588 :       if (ismpzero(gel(x,1))) /* avoid cancellation */
     624         231 :         return gerepilecopy(av, mulcxI(gasin(gel(x,2), prec)));
     625       14357 :       d = gsqrt(gaddsg(1,gsqr(x)), prec); /* Re(d) >= 0 */
     626       14357 :       a = gadd(d, x);
     627       14357 :       b = gsub(d, x);
     628             :       /* avoid cancellation as much as possible */
     629       14357 :       if (gprecision(a) < gprecision(b))
     630           7 :         y = gneg(glog(b,prec));
     631             :       else
     632       14350 :         y = glog(a,prec);
     633       14357 :       return gerepileupto(av, y); /* log (x + sqrt(1+x^2)) */
     634             :     }
     635         196 :     default:
     636         196 :       av = avma; if (!(y = toser_i(x))) break;
     637         168 :       if (gequal0(y)) return gerepilecopy(av, y);
     638         161 :       if (valp(y) < 0) pari_err_DOMAIN("asinh","valuation", "<", gen_0, x);
     639         154 :       p1 = gaddsg(1,gsqr(y));
     640         154 :       if (gequal0(p1))
     641             :       {
     642          14 :         GEN t = PiI2n(-1,prec);
     643          14 :         if ( gsigne(imag_i(gel(y,2))) < 0 ) setsigne(gel(t,2), -1);
     644          14 :         return gerepileupto(av, scalarser(t, varn(y), valp(p1)>>1));
     645             :       }
     646         140 :       p1 = gdiv(derivser(y), gsqrt(p1,prec));
     647         140 :       a = integser(p1);
     648         140 :       if (!valp(y)) a = gadd(a, gasinh(gel(y,2),prec));
     649         140 :       return gerepileupto(av, a);
     650             :   }
     651          28 :   return trans_eval("asinh",gasinh,x,prec);
     652             : }
     653             : /********************************************************************/
     654             : /**                                                                **/
     655             : /**                     AREA HYPERBOLIC COSINE                     **/
     656             : /**                                                                **/
     657             : /********************************************************************/
     658             : 
     659             : /* |x| >= 1, return ach(|x|) */
     660             : static GEN
     661         742 : mpacosh(GEN x)
     662             : {
     663         742 :   long lx = realprec(x), e;
     664         742 :   GEN z, res = cgetr(lx);
     665         742 :   pari_sp av = avma;
     666         742 :   e = expo(signe(x) > 0? subrs(x,1): addrs(x,1));
     667         742 :   if (e == -(long)HIGHEXPOBIT)
     668           0 :     return gc_const((pari_sp)(res + lx), real_0_bit(- bit_prec(x) >> 1));
     669         742 :   if (e < -5) x = rtor(x, realprec(x) + nbits2extraprec(-e));
     670         742 :   z = logr_abs( addrr_sign(x, 1, sqrtr( subrs(sqrr(x), 1) ), 1) );
     671         742 :   affrr(z, res); return gc_const(av, res);
     672             : }
     673             : 
     674             : GEN
     675        8001 : gacosh(GEN x, long prec)
     676             : {
     677             :   pari_sp av;
     678             :   GEN y;
     679             : 
     680        8001 :   switch(typ(x))
     681             :   {
     682         280 :     case t_REAL: {
     683         280 :       long s = signe(x), e = expo(x);
     684             :       GEN a, b;
     685         280 :       if (s > 0 && e >= 0) return mpacosh(x);
     686             :       /* x < 1 */
     687         147 :       y = cgetg(3,t_COMPLEX); a = gen_0;
     688         147 :       if (s == 0) b = acos0(e);
     689         140 :       else if (e < 0) b = mpacos(x); /* -1 < x < 1 */
     690             :       else {
     691          91 :         if (!absrnz_equal1(x)) a = mpacosh(x);
     692          91 :         b = mppi(realprec(x));
     693             :       }
     694         147 :       gel(y,1) = a;
     695         147 :       gel(y,2) = b; return y;
     696             :     }
     697        7644 :     case t_COMPLEX: {
     698             :       GEN a, b, d;
     699        7644 :       if (ismpzero(gel(x,2))) return gacosh(gel(x,1), prec);
     700        7413 :       av = avma;
     701        7413 :       d = gsqrt(gaddsg(-1,gsqr(x)), prec); /* Re(d) >= 0 */
     702        7413 :       a = gadd(x, d);
     703        7413 :       b = gsub(x, d);
     704             :       /* avoid cancellation as much as possible */
     705        7413 :       if (gprecision(a) < gprecision(b))
     706          14 :         y = glog(b,prec);
     707             :       else
     708        7399 :         y = glog(a,prec);
     709             :       /* y = \pm log(x + sqrt(x^2-1)) */
     710        7413 :       if (gsigne(real_i(y)) < 0) y = gneg(y);
     711        7413 :       return gerepileupto(av, y);
     712             :     }
     713          77 :     default: {
     714             :       GEN a, d;
     715             :       long v;
     716          77 :       av = avma; if (!(y = toser_i(x))) break;
     717          49 :       v = valp(y);
     718          49 :       if (v < 0) pari_err_DOMAIN("acosh","valuation", "<", gen_0, x);
     719          42 :       if (gequal0(y))
     720             :       {
     721           7 :         if (!v) return gerepilecopy(av, y);
     722           7 :         return gerepileupto(av, gadd(y, PiI2n(-1, prec)));
     723             :       }
     724          35 :       d = gsubgs(gsqr(y),1);
     725          35 :       if (gequal0(d)) { set_avma(av); return zeroser(varn(y), valp(d)>>1); }
     726          21 :       d = gdiv(derivser(y), gsqrt(d,prec));
     727          21 :       a = integser(d);
     728          21 :       if (v)
     729           7 :         d = PiI2n(-1, prec); /* I Pi/2 */
     730             :       else
     731             :       {
     732          14 :         d = gel(y,2); if (gequal1(d)) return gerepileupto(av,a);
     733           7 :         d = gacosh(d, prec);
     734             :       }
     735          14 :       return gerepileupto(av, gadd(d,a));
     736             :     }
     737             :   }
     738          28 :   return trans_eval("acosh",gacosh,x,prec);
     739             : }
     740             : /********************************************************************/
     741             : /**                                                                **/
     742             : /**                     AREA HYPERBOLIC TANGENT                    **/
     743             : /**                                                                **/
     744             : /********************************************************************/
     745             : 
     746             : /* |x| < 1, x != 0 */
     747             : static GEN
     748         126 : mpatanh(GEN x)
     749             : {
     750         126 :   pari_sp av = avma;
     751         126 :   long e, s = signe(x);
     752             :   GEN z;
     753         126 :   z = s > 0? subsr(1,x): addsr(1,x); e = expo(z);
     754         126 :   if (e < -5)
     755             :   {
     756          28 :     x = rtor(x, realprec(x) + nbits2extraprec(-e)-1);
     757          28 :     z = s > 0? subsr(1,x): addsr(1,x); e = expo(z);
     758             :   }
     759         126 :   z = invr(z); shiftr_inplace(z, 1); /* 2/(1-|x|) */
     760         126 :   z = logr_abs( addrs(z,-1) ); if (s < 0) togglesign(z);
     761         126 :   shiftr_inplace(z, -1); return gerepileuptoleaf(av, z);
     762             : }
     763             : 
     764             : GEN
     765       15631 : gatanh(GEN x, long prec)
     766             : {
     767             :   long sx;
     768             :   pari_sp av;
     769             :   GEN a, y, z;
     770             : 
     771       15631 :   switch(typ(x))
     772             :   {
     773         511 :     case t_REAL:
     774         511 :       sx = signe(x);
     775         511 :       if (!sx) return real_0_bit(expo(x));
     776         504 :       if (expo(x) < 0) return mpatanh(x);
     777             : 
     778         378 :       y = cgetg(3,t_COMPLEX);
     779         378 :       av = avma;
     780         378 :       z = subrs(x,1);
     781         378 :       if (!signe(z)) pari_err_DOMAIN("atanh", "argument", "=", gen_1, x);
     782         364 :       z = invr(z); shiftr_inplace(z, 1); /* 2/(x-1)*/
     783         364 :       z = addrs(z,1);
     784         364 :       if (!signe(z)) pari_err_DOMAIN("atanh", "argument", "=", gen_m1, x);
     785         350 :       z = logr_abs(z);
     786         350 :       shiftr_inplace(z, -1); /* (1/2)log((1+x)/(x-1)) */
     787         350 :       gel(y,1) = gerepileuptoleaf(av, z);
     788         350 :       gel(y,2) = Pi2n(-1, realprec(x));
     789         350 :       if (sx > 0) togglesign(gel(y,2));
     790         350 :       return y;
     791             : 
     792       15071 :     case t_COMPLEX: /* 2/(1-z) - 1 = (1+z) / (1-z) */
     793       15071 :       if (ismpzero(gel(x,2))) return gatanh(gel(x,1), prec);
     794       14616 :       av = avma; z = glog( gaddgs(gdivsg(2,gsubsg(1,x)),-1), prec );
     795       14616 :       return gerepileupto(av, gmul2n(z,-1));
     796             : 
     797          49 :     default:
     798          49 :       av = avma; if (!(y = toser_i(x))) break;
     799          28 :       if (valp(y) < 0) pari_err_DOMAIN("atanh","valuation", "<", gen_0, x);
     800          21 :       z = gdiv(derivser(y), gsubsg(1,gsqr(y)));
     801          14 :       a = integser(z);
     802          14 :       if (!valp(y)) a = gadd(a, gatanh(gel(y,2),prec));
     803          14 :       return gerepileupto(av, a);
     804             :   }
     805          21 :   return trans_eval("atanh",gatanh,x,prec);
     806             : }
     807             : /********************************************************************/
     808             : /**                                                                **/
     809             : /**                         EULER'S GAMMA                          **/
     810             : /**                                                                **/
     811             : /********************************************************************/
     812             : /* 0 < a < b */
     813             : static GEN
     814       25780 : mulu_interval_step_i(ulong a, ulong b, ulong step)
     815             : {
     816             :   ulong k, l, N, n;
     817             :   long lx;
     818             :   GEN x;
     819             : 
     820       25780 :   n = 1 + (b-a) / step;
     821       25780 :   b -= (b-a) % step;
     822             :   /* step | b-a */
     823       25780 :   lx = 1; x = cgetg(2 + n/2, t_VEC);
     824       25780 :   N = b + a;
     825       25780 :   for (k = a;; k += step)
     826             :   {
     827      128138 :     l = N - k; if (l <= k) break;
     828      102358 :     gel(x,lx++) = muluu(k,l);
     829             :   }
     830       25780 :   if (l == k) gel(x,lx++) = utoipos(k);
     831       25780 :   setlg(x, lx); return x;
     832             : }
     833             : static GEN
     834       99605 : _mul(void *data, GEN x, GEN y)
     835             : {
     836       99605 :   long prec = (long)data;
     837             :   /* switch to t_REAL ? */
     838       99605 :   if (typ(x) == t_INT && lgefint(x) > prec) x = itor(x, prec);
     839       99605 :   if (typ(y) == t_INT && lgefint(y) > prec) y = itor(y, prec);
     840       99605 :   return mpmul(x, y);
     841             : }
     842             : static GEN
     843       25780 : mulu_interval_step_prec(long l, long m, long s, long prec)
     844             : {
     845       25780 :   GEN v = mulu_interval_step_i(l, m, s);
     846       25780 :   return gen_product(v, (void*)prec, &_mul);
     847             : }
     848             : 
     849             : /* x * (i*(i+1)) */
     850             : static GEN
     851     8487666 : muliunu(GEN x, ulong i)
     852             : {
     853     8487666 :   if (i & HIGHMASK) /* i(i+1) >= 2^BITS_IN_LONG*/
     854           0 :     return muliu(muliu(x, i), i+1);
     855             :   else
     856     8487666 :     return muliu(x, i*(i+1));
     857             : }
     858             : /* x / (i*(i+1)) */
     859             : GEN
     860    25219072 : divrunu(GEN x, ulong i)
     861             : {
     862    25219072 :   if (i & HIGHMASK) /* i(i+1) >= 2^BITS_IN_LONG*/
     863           0 :     return divru(divru(x, i), i+1);
     864             :   else
     865    25219072 :     return divru(x, i*(i+1));
     866             : }
     867             : /* x / (i*(i+1)) */
     868             : GEN
     869      779559 : divgunu(GEN x, ulong i)
     870             : {
     871             : #ifdef LONG_IS_64BIT
     872      670128 :   if (i < 3037000500L) /* i(i+1) < 2^63 */
     873             : #else
     874      109431 :   if (i < 46341L) /* i(i+1) < 2^31 */
     875             : #endif
     876      779559 :     return gdivgs(x, i*(i+1));
     877             :   else
     878           0 :     return gdivgs(gdivgs(x, i), i+1);
     879             : }
     880             : 
     881             : /* arg(s+it) */
     882             : double
     883      236970 : darg(double s, double t)
     884             : {
     885             :   double x;
     886      236970 :   if (!t) return (s>0)? 0.: M_PI;
     887      206129 :   if (!s) return (t>0)? M_PI/2: -M_PI/2;
     888      206122 :   x = atan(t/s);
     889             :   return (s>0)? x
     890      206122 :               : ((t>0)? x+M_PI : x-M_PI);
     891             : }
     892             : 
     893             : void
     894      236968 : dcxlog(double s, double t, double *a, double *b)
     895             : {
     896      236968 :   *a = log(s*s + t*t) / 2; /* log |s| = Re(log(s)) */
     897      236968 :   *b = darg(s,t);          /* Im(log(s)) */
     898      236970 : }
     899             : 
     900             : double
     901       16212 : dabs(double s, double t) { return sqrt( s*s + t*t ); }
     902             : double
     903        4102 : dnorm(double s, double t) { return s*s + t*t; }
     904             : 
     905             : #if 0
     906             : /* x, z t_REAL. Compute unique x in ]-z,z] congruent to x mod 2z */
     907             : static GEN
     908             : red_mod_2z(GEN x, GEN z)
     909             : {
     910             :   GEN Z = gmul2n(z, 1), d = subrr(z, x);
     911             :   /* require little accuracy */
     912             :   if (!signe(d)) return x;
     913             :   setprec(d, nbits2prec(expo(d) - expo(Z)));
     914             :   return addrr(mulir(floorr(divrr(d, Z)), Z), x);
     915             : }
     916             : #endif
     917             : 
     918             : static GEN
     919        6167 : negeuler(long prec) { GEN g = mpeuler(prec); setsigne(g, -1); return g; }
     920             : /* lngamma(1+z) = -Euler*z + sum_{i > 1} zeta(i)/i (-z)^i
     921             :  * at relative precision prec, |z| < 1 is small */
     922             : static GEN
     923        2475 : lngamma1(GEN z, long prec)
     924             : { /* sum_{i > l} |z|^(i-1) = |z|^l / (1-|z|) < 2^-B
     925             :    * for l > (B+1) / |log2(|z|)| */
     926        2475 :   long i, l = ceil((bit_accuracy(prec) + 1) / - dbllog2(z));
     927             :   GEN s, vz;
     928             : 
     929        2475 :   if (l <= 1) return gmul(negeuler(prec), z);
     930        2300 :   vz = constzeta(l, prec);
     931       18910 :   for (i = l, s = gen_0; i > 0; i--)
     932             :   {
     933       16610 :     GEN c = divru(gel(vz,i), i);
     934       16610 :     if (odd(i)) setsigne(c, -1);
     935       16610 :     s = gadd(gmul(s,z), c);
     936             :   }
     937        2300 :   return gmul(z, s);
     938             : }
     939             : /* B_i / (i(i-1)), i even. Sometimes NOT reduced (but gadd/gmul won't care)!*/
     940             : static GEN
     941     8487511 : bern_unu(long i)
     942     8487511 : { GEN B = bernfrac(i); return mkfrac(gel(B,1), muliunu(gel(B,2), i-1)); }
     943             : /* B_i / i, i even. Sometimes NOT reduced (but gadd/gmul won't care)!*/
     944             : static GEN
     945      206186 : bern_u(long i)
     946      206186 : { GEN B = bernfrac(i); return mkfrac(gel(B,1), muliu(gel(B,2), i)); }
     947             : /* sum_{i > 0} B_{2i}/(2i(2i-1)) * a^(i-1) */
     948             : static GEN
     949      229154 : lngamma_sum(GEN a, long N)
     950             : {
     951      229154 :   pari_sp av = avma;
     952      229154 :   GEN S = bern_unu(2*N);
     953             :   long i;
     954     8488331 :   for (i = 2*N-2; i > 0; i -= 2)
     955             :   {
     956     8259247 :     S = gadd(bern_unu(i), gmul(a,S));
     957     8259279 :     if (gc_needed(av,3))
     958             :     {
     959           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"gamma: i = %ld", i);
     960           0 :       S = gerepileupto(av, S);
     961             :     }
     962             :   }
     963      229084 :   return S;
     964             : }
     965             : /* sum_{i > 0} B_{2i}/(2i) * a^i */
     966             : static GEN
     967        4109 : psi_sum(GEN a, long N)
     968             : {
     969        4109 :   pari_sp av = avma;
     970        4109 :   GEN S = bern_u(2*N);
     971             :   long i;
     972      206186 :   for (i = 2*N-2; i > 0; i -= 2)
     973             :   {
     974      202077 :     S = gadd(bern_u(i), gmul(a,S));
     975      202077 :     if (gc_needed(av,3))
     976             :     {
     977           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"psi: i = %ld", i);
     978           0 :       S = gerepileupto(av, S);
     979             :     }
     980             :   }
     981        4109 :   return gmul(a,S);
     982             : }
     983             : static GEN
     984      231966 : cxgamma(GEN s0, int dolog, long prec)
     985             : {
     986             :   GEN s, a, y, res, sig, tau, p1, nnx, pi, pi2, sqrtpi2;
     987             :   long i, lim, N, esig, et;
     988             :   pari_sp av, av2;
     989      231966 :   int funeq = 0;
     990             :   pari_timer T;
     991             : 
     992      231966 :   if (DEBUGLEVEL>5) timer_start(&T);
     993      231966 :   s = trans_fix_arg(&prec,&s0,&sig,&tau,&av,&res);
     994             : 
     995      231955 :   esig = expo(sig);
     996      231955 :   et = signe(tau)? expo(tau): 0;
     997      231955 :   if ((signe(sig) <= 0 || esig < -1) && et <= 16)
     998             :   { /* s <--> 1-s */
     999       30737 :     funeq = 1; s = gsubsg(1, s); sig = real_i(s);
    1000             :   }
    1001             : 
    1002             :   /* find "optimal" parameters [lim, N] */
    1003      231955 :   if (esig > 300 || et > 300)
    1004          35 :   { /* |s| is HUGE ! Play safe and avoid inf / NaN */
    1005             :     GEN S, iS, l2, la, u;
    1006             :     double logla, l;
    1007             : 
    1008          36 :     S = gprec_w(s,LOWDEFAULTPREC);
    1009             :     /* l2 ~ |lngamma(s))|^2 */
    1010          35 :     l2 = gnorm(gmul(S, glog(S, LOWDEFAULTPREC)));
    1011          35 :     l = (prec2nbits_mul(prec, M_LN2) - rtodbl(glog(l2,LOWDEFAULTPREC))/2) / 2.;
    1012          35 :     if (l < 0) l = 0.;
    1013             : 
    1014          35 :     iS = imag_i(S);
    1015          35 :     if (et > 0 && l > 0)
    1016          21 :     {
    1017          21 :       GEN t = gmul(iS, dbltor(M_PI / l)), logt = glog(t,LOWDEFAULTPREC);
    1018          21 :       la = gmul(t, logt);
    1019          21 :       if      (gcmpgs(la, 3) < 0)   { logla = log(3.); la = stoi(3); }
    1020          14 :       else if (gcmpgs(la, 150) > 0) { logla = rtodbl(logt); la = t; }
    1021           7 :       else logla = rtodbl(mplog(la));
    1022             :     }
    1023             :     else
    1024             :     {
    1025          14 :       logla = log(3.); la = stoi(3);
    1026             :     }
    1027          35 :     lim = (long)ceil(l / (1.+ logla));
    1028          35 :     if (lim == 0) lim = 1;
    1029             : 
    1030          35 :     u = gmul(la, dbltor((lim-0.5)/M_PI));
    1031          35 :     l2 = gsub(gsqr(u), gsqr(iS));
    1032          35 :     if (signe(l2) > 0)
    1033             :     {
    1034          14 :       l2 = gsub(gsqrt(l2,3), sig);
    1035          14 :       if (signe(l2) > 0) N = itos( gceil(l2) ); else N = 1;
    1036             :     }
    1037             :     else
    1038          21 :       N = 1;
    1039             :   }
    1040             :   else
    1041             :   { /* |s| is moderate. Use floats  */
    1042      231919 :     double ssig = rtodbl(sig);
    1043      231918 :     double st = typ(s) == t_REAL? 0.0: rtodbl(imag_i(s));
    1044             :     double la, l,l2,u,v, rlogs, ilogs;
    1045             : 
    1046      231925 :     if (fabs(ssig-1) + fabs(st) < 1e-16)
    1047             :     { /* s ~ 1: loggamma(1+u) ~ - Euler * u, cancellation */
    1048        2807 :       if (funeq) /* s0 ~ 0: use lngamma(s0)+log(s0) = lngamma(s0+1) */
    1049             :       {
    1050          56 :         if (dolog)
    1051           0 :           y = gsub(lngamma1(s0,prec), glog(s0,prec));
    1052             :         else
    1053          56 :           y = gdiv(gexp(lngamma1(s0,prec), prec), s0);
    1054             :       }
    1055             :       else
    1056             :       {
    1057        3847 :         if (isint1(s0)) { set_avma(av); return dolog? real_0(prec): real_1(prec); }
    1058        1040 :         y = lngamma1(gsubgs(s0,1),prec);
    1059        1040 :         if (!dolog) y = gexp(y,prec);
    1060             :       }
    1061        1096 :       set_avma(av); return affc_fixlg(y, res);
    1062             :     }
    1063      229118 :     dcxlog(ssig,st, &rlogs,&ilogs);
    1064             :     /* Re (s - 1/2) log(s) */
    1065      229121 :     u = (ssig - 0.5)*rlogs - st * ilogs;
    1066             :     /* Im (s - 1/2) log(s) */
    1067      229121 :     v = (ssig - 0.5)*ilogs + st * rlogs;
    1068             :     /* l2 = | (s - 1/2) log(s) - s + log(2Pi)/2 |^2 ~ |lngamma(s))|^2 */
    1069      229121 :     u = u - ssig + log(2.*M_PI)/2;
    1070      229121 :     v = v - st;
    1071      229121 :     l2 = u*u + v*v;
    1072      229121 :     if (l2 < 0.000001) l2 = 0.000001;
    1073      229121 :     l = (prec2nbits_mul(prec, M_LN2) - log(l2)/2) / 2.;
    1074      229122 :     if (l < 0) l = 0.;
    1075             : 
    1076      229122 :     if (st > 1 && l > 0)
    1077       63973 :     {
    1078       63973 :       double t = st * M_PI / l;
    1079       63973 :       la = t * log(t);
    1080       63973 :       if (la < 4.) la = 4.;
    1081       63973 :       if (la > 150) la = t;
    1082             :     }
    1083             :     else
    1084      165149 :       la = 4.; /* heuristic */
    1085      229122 :     lim = (long)ceil(l / (1.+ log(la)));
    1086      229122 :     if (lim == 0) lim = 1;
    1087             : 
    1088      229122 :     u = (lim-0.5) * la / M_PI;
    1089      229122 :     l2 = u*u - st*st;
    1090      229122 :     if (l2 > 0)
    1091             :     {
    1092      219557 :       double t = ceil(sqrt(l2) - ssig);
    1093      219557 :       N = (t < 1)? 1: (long)t;
    1094      219557 :       if (N < 1) N = 1;
    1095             :     }
    1096             :     else
    1097        9565 :       N = 1;
    1098             :   }
    1099      229157 :   if (DEBUGLEVEL>5) err_printf("lim, N: [%ld, %ld]\n",lim,N);
    1100      229163 :   incrprec(prec);
    1101             : 
    1102      229163 :   av2 = avma;
    1103      229163 :   y = s;
    1104      229163 :   if (typ(s0) == t_INT)
    1105             :   {
    1106        2381 :     ulong ss = itou_or_0(s0);
    1107        2381 :     if (signe(s0) <= 0)
    1108           0 :       pari_err_DOMAIN("gamma","argument", "=",
    1109             :                        strtoGENstr("non-positive integer"), s0);
    1110        2388 :     if (!ss || ss + (ulong)N < ss) {
    1111           7 :       for (i=1; i < N; i++)
    1112             :       {
    1113           0 :         y = mulri(y, addiu(s0, i));
    1114           0 :         if (gc_needed(av2,3))
    1115             :         {
    1116           0 :           if(DEBUGMEM>1) pari_warn(warnmem,"gamma");
    1117           0 :           y = gerepileuptoleaf(av2, y);
    1118             :         }
    1119             :       }
    1120             :     } else {
    1121       31550 :       for (i=1; i < N; i++)
    1122             :       {
    1123       29176 :         y = mulru(y, ss + i);
    1124       29176 :         if (gc_needed(av2,3))
    1125             :         {
    1126           0 :           if(DEBUGMEM>1) pari_warn(warnmem,"gamma");
    1127           0 :           y = gerepileuptoleaf(av2, y);
    1128             :         }
    1129             :       }
    1130             :     }
    1131        2381 :     if (dolog) y = logr_abs(y);
    1132             :   }
    1133             :   else
    1134             :   { /* Compute lngamma mod 2 I Pi */
    1135      226782 :     GEN sq = gsqr(s);
    1136      226779 :     pari_sp av3 = avma;
    1137     4765767 :     for (i = 1; i < N - 1; i += 2)
    1138             :     {
    1139     4538991 :       y = gmul(y, gaddsg(i*(i + 1), gadd(gmulsg(2*i + 1, s), sq)));
    1140     4538985 :       if (gc_needed(av2,3))
    1141             :       {
    1142           0 :         if(DEBUGMEM>1) pari_warn(warnmem,"gamma");
    1143           0 :         y = gerepileupto(av3, y);
    1144             :       }
    1145             :     }
    1146      226776 :     if (!odd(N)) y = gmul(y, gaddsg(N - 1, s));
    1147      226776 :     if (dolog)
    1148             :     {
    1149       18332 :       if (typ(s) == t_REAL) y = logr_abs(y);
    1150             :       else
    1151             :       { /* fix imaginary part */
    1152         273 :         long prec0 = LOWDEFAULTPREC;
    1153         273 :         GEN s0 = gprec_w(s, prec0), y0 = s0, k;
    1154         273 :         y0 = garg(y0, prec0); /* Im log(s) at low accuracy */
    1155       10742 :         for (i=1; i < N; i++) y0 = gadd(y0, garg(gaddgs(s0,i), prec0));
    1156         273 :         y = glog(y, prec);
    1157         273 :         k = ground( gdiv(gsub(y0, imag_i(y)), Pi2n(1,prec0)) );
    1158         273 :         if (signe(k)) y = gadd(y, mulcxI(mulir(k, Pi2n(1, prec))));
    1159             :       }
    1160             :     }
    1161             :   }
    1162      229157 :   if (DEBUGLEVEL>5) timer_printf(&T,"product from 0 to N-1");
    1163      229157 :   constbern(lim);
    1164      229158 :   nnx = gaddgs(s, N); a = ginv(nnx);
    1165      229152 :   p1 = gsub(gmul(gsub(nnx, ghalf), glog(nnx,prec)), nnx);
    1166      229147 :   p1 = gadd(p1, gmul(a, lngamma_sum(gsqr(a), lim)));
    1167      229152 :   if (DEBUGLEVEL>5) timer_printf(&T,"Bernoulli sum");
    1168             : 
    1169      229152 :   pi = mppi(prec); pi2 = shiftr(pi, 1); sqrtpi2 = sqrtr(pi2);
    1170             : 
    1171      229148 :   if (dolog)
    1172             :   {
    1173       18578 :     if (funeq)
    1174             :     { /* recall that s = 1 - s0 */
    1175          14 :       GEN T = shiftr(sqrtpi2,-1); /* sqrt(2Pi)/2 */
    1176          14 :       if (typ(s) != t_REAL)
    1177             :       {
    1178             :         /* We compute log(sin(Pi s0)) so that it has branch cuts along
    1179             :         * (-oo, 0] and [1, oo).  To do this in a numerically stable way
    1180             :         * we must compute the log first then mangle its imaginary part.
    1181             :         * The rounding operation below is stable because we're rounding
    1182             :         * a number which is already within 1/4 of an integer. */
    1183             : 
    1184             :         /* z = log( sin(Pi s0) / (sqrt(2Pi)/2) ) */
    1185           7 :         GEN z = glog(gdiv(gsin(gmul(pi,s0),prec), T), prec);
    1186             :         /* b = (2 Re(s) - 1) / 4 */
    1187           7 :         GEN b = shiftr(subrs(shiftr(sig, 1), 1), -2);
    1188           7 :         y = gsub(y, z);
    1189           7 :         if (gsigne(imag_i(s)) > 0) togglesign(b);
    1190             :         /* z = 2Pi round( Im(z)/2Pi - b ) */
    1191           7 :         z = gmul(roundr(gsub(gdiv(imag_i(z), pi2), b)), pi2);
    1192           7 :         if (signe(z)) { /* y += I*z, z a t_REAL */
    1193           0 :           if (typ(y) == t_COMPLEX)
    1194           0 :             gel(y,2) = gadd(gel(y,2), z);
    1195             :           else
    1196           0 :             y = mkcomplex(y, z);
    1197             :         }
    1198             :       }
    1199             :       else
    1200             :       { /* s0 < 0, formula simplifies: imag(lngamma(s0)) = - Pi * floor(s0) */
    1201           7 :         GEN z = logr_abs(divrr(mpsin(gmul(pi,s0)), T));
    1202           7 :         y = gsub(y, z);
    1203           7 :         y = mkcomplex(y, mulri(pi, gfloor(s0)));
    1204             :       }
    1205          14 :       p1 = gneg(p1);
    1206             :     }
    1207             :     else /* y --> sqrt(2Pi) / y */
    1208       18564 :       y = gsub(logr_abs(sqrtpi2), y);
    1209       18578 :     y = gadd(p1, y);
    1210             :   }
    1211             :   else
    1212             :   {
    1213      210570 :     if (funeq)
    1214             :     { /* y --> y Pi/(sin(Pi s) * sqrt(2Pi)) = y sqrt(Pi/2)/sin(Pi s) */
    1215       30667 :       y = gdiv(gmul(shiftr(sqrtpi2,-1),y), gsin(gmul(pi,s0), prec));
    1216             :       /* don't use s above: sin(pi s0) = sin(pi s) and the former is
    1217             :        * more accurate, esp. if s0 ~ 0 */
    1218       30667 :       p1 = gneg(p1);
    1219             :     }
    1220             :     else /* y --> sqrt(2Pi) / y */
    1221      179903 :       y = gdiv(sqrtpi2, y);
    1222      210574 :     y = gmul(gexp(p1, prec), y);
    1223             :   }
    1224      229148 :   set_avma(av); return affc_fixlg(y, res);
    1225             : }
    1226             : 
    1227             : /* Theory says n > C * b^1.5 / log(b). Timings:
    1228             :  * b = 64*[1, 2, 3, 4, 5, 6, 7, 10, 20, 30, 50, 100, 200, 500];
    1229             :  * n = [1450, 1930, 2750, 3400, 4070, 5000, 6000, 8800, 26000, 50000, 130000,
    1230             :  *      380000, 1300000, 6000000]; */
    1231             : static long
    1232       34517 : gamma2_n(long prec)
    1233             : {
    1234       34517 :   long b = bit_accuracy(prec);
    1235       34517 :   if (b <=  64) return 1450;
    1236       33831 :   if (b <= 128) return 1930;
    1237       28189 :   if (b <= 192) return 2750;
    1238       16319 :   if (b <= 256) return 3400;
    1239        6848 :   if (b <= 320) return 4070;
    1240        6511 :   if (b <= 384) return 5000;
    1241        3985 :   if (b <= 448) return 6000;
    1242        3779 :   return 10.0 * b * sqrt(b) / log(b);
    1243             : }
    1244             : 
    1245             : /* m even, Gamma((m+1) / 2) */
    1246             : static GEN
    1247       34517 : gammahs(long m, long prec)
    1248             : {
    1249       34517 :   GEN y = cgetr(prec), z;
    1250       34517 :   pari_sp av = avma;
    1251       34517 :   long ma = labs(m);
    1252             : 
    1253       34517 :   if (ma > gamma2_n(prec))
    1254             :   {
    1255           0 :     z = stor(m + 1, prec); shiftr_inplace(z, -1);
    1256           0 :     affrr(cxgamma(z,0,prec), y);
    1257           0 :     set_avma(av); return y;
    1258             :   }
    1259       34517 :   z = sqrtr( mppi(prec) );
    1260       34517 :   if (m)
    1261             :   {
    1262       21567 :     GEN t = mulu_interval_step_prec(1, ma-1, 2, prec + EXTRAPRECWORD);
    1263       21567 :     if (m >= 0) z = mpmul(z,t);
    1264             :     else
    1265             :     {
    1266         203 :       z = mpdiv(z,t);
    1267         203 :       if ((m&3) == 2) setsigne(z,-1);
    1268             :     }
    1269       21567 :     shiftr_inplace(z, -m/2);
    1270             :   }
    1271       34517 :   affrr(z, y); set_avma(av); return y;
    1272             : }
    1273             : GEN
    1274          28 : ggammah(GEN x, long prec)
    1275             : {
    1276          28 :   switch(typ(x))
    1277             :   {
    1278          21 :     case t_INT:
    1279             :     {
    1280          21 :       long k = itos_or_0(x);
    1281          21 :       if (!k && signe(x)) pari_err_OVERFLOW("gamma");
    1282          21 :       return gammahs(k * 2, prec);
    1283             :     }
    1284           7 :     case t_REAL: case t_COMPLEX: case t_PADIC: case t_SER: {
    1285           7 :       pari_sp av = avma;
    1286           7 :       return gerepileupto(av, ggamma(gadd(x,ghalf), prec));
    1287             :     }
    1288             :   }
    1289           0 :   return trans_eval("gammah",ggammah,x,prec);
    1290             : }
    1291             : 
    1292             : /* find n such that n+v_p(n!)>=k p^2/(p-1)^2 */
    1293             : static long
    1294          28 : nboft(long k, long p)
    1295             : {
    1296          28 :   pari_sp av = avma;
    1297             :   long s, n;
    1298             : 
    1299          28 :   if (k <= 0) return 0;
    1300          28 :   k = itou( gceil(gdiv(mului(k, sqru(p)), sqru(p-1))) );
    1301          28 :   set_avma(av);
    1302         406 :   for (s=0, n=0; n+s < k; n++, s += u_lval(n, p));
    1303          28 :   return n;
    1304             : }
    1305             : 
    1306             : /* Using Dwork's expansion, compute \Gamma(px+1)=-\Gamma(px) with x a unit.
    1307             :  * See p-Adic Gamma Functions and Dwork Cohomology, Maurizio Boyarsky
    1308             :  * Transactions of the AMS, Vol. 257, No. 2. (Feb., 1980), pp. 359-369.
    1309             :  * Inspired by a GP script by Fernando Rodriguez-Villegas */
    1310             : static GEN
    1311          28 : gadw(GEN x, long p)
    1312             : {
    1313          28 :   pari_sp ltop = avma;
    1314          28 :   GEN s, t, u = cgetg(p+1, t_VEC);
    1315          28 :   long j, k, kp, n = nboft(precp(x)+valp(x)+1, p);
    1316             : 
    1317          28 :   t = s = gaddsg(1, zeropadic(gel(x,2), n));
    1318          28 :   gel(u, 1) = s;
    1319          28 :   gel(u, 2) = s;
    1320          42 :   for (j = 2; j < p; ++j)
    1321          14 :     gel(u, j+1) = gdivgs(gel(u, j), j);
    1322         378 :   for (k = 1, kp = p; k < n; ++k, kp += p) /* kp = k*p */
    1323             :   {
    1324             :     GEN c;
    1325         350 :     gel(u, 1) = gdivgs(gadd(gel(u, 1), gel(u, p)), kp);
    1326         812 :     for (j = 1; j < p; ++j)
    1327         462 :       gel(u, j+1) = gdivgs(gadd(gel(u, j), gel(u, j+1)), kp + j);
    1328             : 
    1329         350 :     t = gmul(t, gaddgs(x, k-1));
    1330         350 :     c = leafcopy(gel(u,1)); setvalp(c, valp(c) + k); /* c = u[1] * p^k */
    1331         350 :     s = gadd(s, gmul(c, t));
    1332         350 :     if ((k&0xFL)==0) gerepileall(ltop, 3, &u,&s,&t);
    1333             :   }
    1334          28 :   return gneg(s);
    1335             : }
    1336             : 
    1337             : /*Use Dwork expansion*/
    1338             : /*This is a O(p*e*log(pe)) algorithm, should be used when p small
    1339             :  * If p==2 this is a O(pe) algorithm. */
    1340             : static GEN
    1341          28 : Qp_gamma_Dwork(GEN x, long p)
    1342             : {
    1343          28 :   pari_sp ltop = avma;
    1344          28 :   long k = padic_to_Fl(x, p);
    1345             :   GEN p1;
    1346             :   long j;
    1347          28 :   long px = precp(x);
    1348          28 :   if (p==2 && px)
    1349             :   {
    1350          14 :     x = shallowcopy(x);
    1351          14 :     setprecp(x, px+1);
    1352          14 :     gel(x,3) = shifti(gel(x,3),1);
    1353             :   }
    1354          28 :   if (k)
    1355             :   {
    1356          21 :     GEN x_k = gsubgs(x,k);
    1357          21 :     x = gdivgs(x_k, p);
    1358          21 :     p1 = gadw(x, p); if (!odd(k)) p1 = gneg(p1);
    1359          28 :     for (j = 1; j < k; ++j) p1 = gmul(p1, gaddgs(x_k, j));
    1360             :   }
    1361             :   else
    1362           7 :     p1 = gneg(gadw(gdivgs(x, p), p));
    1363          28 :   return gerepileupto(ltop, p1);
    1364             : }
    1365             : 
    1366             : /* Compute Qp_gamma using the definition. This is a O(x*M(log(pe))) algorithm.
    1367             :  * This should be used if x is very small. */
    1368             : static GEN
    1369          49 : Qp_gamma_Morita(long n, GEN p, long e)
    1370             : {
    1371          49 :   pari_sp ltop=avma;
    1372          49 :   GEN p2 = gaddsg((n&1)?-1:1, zeropadic(p, e));
    1373             :   long i;
    1374          49 :   long pp=is_bigint(p)? 0: itos(p);
    1375         154 :   for (i = 2; i < n; i++)
    1376         105 :     if (!pp || i%pp)
    1377             :     {
    1378          63 :       p2 = gmulgs(p2, i);
    1379          63 :       if ((i&0xFL) == 0xFL)
    1380           0 :         p2 = gerepileupto(ltop, p2);
    1381             :     }
    1382          49 :   return gerepileupto(ltop, p2);
    1383             : }
    1384             : 
    1385             : /* x\in\N: Gamma(-x)=(-1)^(1+x+x\p)*Gamma(1+x) */
    1386             : static GEN
    1387          28 : Qp_gamma_neg_Morita(long n, GEN p, long e)
    1388             : {
    1389          28 :   GEN g = ginv(Qp_gamma_Morita(n+1, p, e));
    1390          28 :   return ((n^sdivsi(n,p)) & 1)? g: gneg(g);
    1391             : }
    1392             : 
    1393             : /* p-adic Gamma function for x a p-adic integer */
    1394             : /* If n < p*e : use Morita's definition.
    1395             :  * Else : use Dwork's expansion.
    1396             :  * If both n and p are big : itos(p) will fail.
    1397             :  * TODO: handle p=2 better (Qp_gamma_Dwork is slow for p=2). */
    1398             : GEN
    1399          77 : Qp_gamma(GEN x)
    1400             : {
    1401          77 :   GEN n, m, N, p = gel(x,2);
    1402          77 :   long s, e = precp(x);
    1403          77 :   if (absequaliu(p, 2) && e == 2) e = 1;
    1404          77 :   if (valp(x) < 0) pari_err_DOMAIN("gamma","v_p(x)", "<", gen_0, x);
    1405          77 :   n = gtrunc(x);
    1406          77 :   m = gtrunc(gneg(x));
    1407          77 :   N = cmpii(n,m)<=0?n:m;
    1408          77 :   s = itos_or_0(N);
    1409          77 :   if (s && cmpsi(s, muliu(p,e)) < 0) /* s < p*e */
    1410          49 :     return (N == n) ? Qp_gamma_Morita(s,p,e): Qp_gamma_neg_Morita(s,p,e);
    1411          28 :   return Qp_gamma_Dwork(x, itos(p));
    1412             : }
    1413             : 
    1414             : /* gamma(1+x) - 1, |x| < 1 is "small" */
    1415             : GEN
    1416        1211 : ggamma1m1(GEN x, long prec) { return gexpm1(lngamma1(x, prec), prec); }
    1417             : 
    1418             : /* lngamma(y) with 0 constant term, using (lngamma y)' = y' psi(y) */
    1419             : static GEN
    1420       19040 : serlngamma0(GEN y, long prec)
    1421             : {
    1422             :   GEN t;
    1423       19040 :   if (valp(y)) pari_err_DOMAIN("lngamma","valuation", "!=", gen_0, y);
    1424       19033 :   t = derivser(y);
    1425             :   /* don't compute psi if y'=0 */
    1426       19033 :   if (signe(t)) t = gmul(t, gpsi(y,prec));
    1427       19033 :   return integser(t);
    1428             : }
    1429             : 
    1430             : static GEN
    1431       19005 : sergamma(GEN y, long prec)
    1432             : {
    1433             :   GEN z, y0, Y;
    1434       19005 :   if (lg(y) == 2) pari_err_DOMAIN("gamma", "argument", "=", gen_0,y);
    1435             :   /* exp(lngamma) */
    1436       18998 :   if (valp(y) > 0) return gdiv(gexp(glngamma(gaddgs(y,1),prec),prec),y);
    1437       18746 :   y0 = simplify_shallow(gel(y,2));
    1438       18746 :   z = NULL; Y = y;
    1439       18746 :   if (isint(y0, &y0))
    1440             :   { /* fun eq. avoids log singularity of lngamma at negative ints */
    1441        9632 :     long s = signe(y0);
    1442             :     /* possible if y[2] is an inexact 0 */
    1443        9632 :     if (!s) return gdiv(gexp(glngamma(gaddgs(y,1),prec),prec),y);
    1444        9625 :     if (signe(y0) < 0) { Y = gsubsg(1, y); y0 = subsi(1, y0); }
    1445        9625 :     if (abscmpiu(y0, 50) < 0) z = mpfact(itos(y0)-1); /* more precise */
    1446             :   }
    1447       18739 :   if (!z) z = ggamma(y0,prec);
    1448       18739 :   z = gmul(z, gexp(serlngamma0(Y,prec),prec));
    1449       18739 :   if (Y != y)
    1450             :   {
    1451          77 :     GEN pi = mppi(prec);
    1452          77 :     z = gdiv(mpodd(y0)? pi: negr(pi),
    1453             :              gmul(z, gsin(gmul(pi,serchop0(y)), prec)));
    1454             :   }
    1455       18739 :   return z;
    1456             : }
    1457             : 
    1458             : static GEN
    1459        8568 : sqrtu(ulong a, long prec) { return sqrtr_abs(utor(a, prec)); }
    1460             : static GEN
    1461         112 : cbrtu(ulong a, long prec) { return sqrtnr_abs(utor(a, prec), 3); }
    1462             : /* N | 6 */
    1463             : static GEN
    1464        5453 : ellkprime(long N, GEN s2, GEN s3)
    1465             : {
    1466             :   GEN z;
    1467        5453 :   switch(N)
    1468             :   {
    1469        1918 :     case 1: return shiftr(s2, -1);
    1470          42 :     case 2: return sqrtr_abs(shiftr(subrs(s2,1), 1));
    1471        3395 :     case 3: return shiftr(mulrr(s2, addrs(s3, 1)), -2);
    1472          98 :     default: /* 6 */
    1473          98 :       z = mulrr(subrr(s3,s2), subsr(2,s3));
    1474          98 :       return mulrr(addsr(2,s2), sqrtr_abs(z));
    1475             :   }
    1476             : }
    1477             : 
    1478             : static GEN
    1479        5453 : ellKk(long N, GEN s2, GEN s3, long prec)
    1480        5453 : { return gdiv(Pi2n(-1,prec), agm(ellkprime(N,s2,s3), gen_1, prec)); }
    1481             : 
    1482             : /* Gamma(1/3) */
    1483             : static GEN
    1484        3297 : G3(GEN s2, GEN s3, long prec)
    1485             : {
    1486        3297 :   GEN A = ellKk(3, s2,s3, prec), pi = mppi(prec);
    1487        3297 :   A = shiftr(divrs(powrs(mulrr(pi, A), 12), 27), 28);
    1488        3297 :   return sqrtnr_abs(A, 36);
    1489             : }
    1490             : /* Gamma(1/4) */
    1491             : static GEN
    1492        1778 : G4(GEN s2, long prec)
    1493             : {
    1494        1778 :   GEN A = ellKk(1, s2,NULL, prec), pi = mppi(prec);
    1495        1778 :   return shiftr(sqrtr_abs(mulrr(sqrtr_abs(pi), A)), 1);
    1496             : }
    1497             : 
    1498             : /* Gamma(n / 24), n = 1,5,7,11 */
    1499             : static GEN
    1500          98 : Gn24(long n, GEN s2, GEN s3, long prec)
    1501             : {
    1502          98 :   GEN A, B, C, t, t1, t2, t3, t4, pi = mppi(prec);
    1503          98 :   A = ellKk(1, s2,s3, prec);
    1504          98 :   B = ellKk(3, s2,s3, prec);
    1505          98 :   C = ellKk(6, s2,s3, prec);
    1506          98 :   t1 = sqrtr_abs(mulur(3, addsr(2, s3)));
    1507          98 :   t2 = sqrtr_abs(divrr(s3, pi));
    1508          98 :   t2 = mulrr(t2, shiftr(mulrr(addrr(s2,s3), A), 2));
    1509          98 :   t3 = mulrr(divur(3,pi), sqrr(B));
    1510          98 :   t3 = mulrr(addsr(2,s2), sqrtnr_abs(shiftr(powrs(t3, 3), 8), 9));
    1511          98 :   t4 = mulrr(mulrr(addsr(1, s2), subrr(s3, s2)), subsr(2, s3));
    1512          98 :   t4 = mulrr(mulrr(mulur(384, t4), pi), sqrr(C));
    1513          98 :   switch (n)
    1514             :   {
    1515          56 :     case 1: t = mulrr(mulrr(t1, t2), mulrr(t3, t4)); break;
    1516          14 :     case 5: t = divrr(mulrr(t2, t4), mulrr(t1, t3)); break;
    1517          14 :     case 7: t = divrr(mulrr(t3, t4), mulrr(t1, t2)); break;
    1518          14 :     default:t = divrr(mulrr(t1, t4), mulrr(t2, t3)); break;
    1519             :   }
    1520          98 :   return sqrtnr_abs(t, 4);
    1521             : }
    1522             : /* sin(x/2) = sqrt((1-c) / 2) > 0 given c = cos(x) */
    1523             : static GEN
    1524          28 : sinx2(GEN c)
    1525          28 : { c = subsr(1, c); shiftr_inplace(c,-1); return sqrtr_abs(c); }
    1526             : /* sin(Pi/12), given sqrt(3) */
    1527             : static GEN
    1528          21 : sin12(GEN s3)
    1529          21 : { GEN t = subsr(2, s3); shiftr_inplace(t, -2); return sqrtr_abs(t); }
    1530             : /* cos(Pi/12) = sin(5Pi/12), given sqrt(3) */
    1531             : static GEN
    1532          21 : cos12(GEN s3)
    1533          21 : { GEN t = addsr(2, s3); shiftr_inplace(t, -2); return sqrtr_abs(t); }
    1534             : /* 0 < n < d, (n, d) = 1, 2 < d | 24 */
    1535             : static GEN
    1536        5173 : gammafrac24_s(long n, long d, long prec)
    1537             : {
    1538        5173 :   GEN A, B, s2, s3, pi = mppi(prec);
    1539        5173 :   s2 = sqrtu(2, prec);
    1540        5173 :   s3 = d % 3? NULL: sqrtu(3, prec);
    1541        5173 :   switch(d)
    1542             :   {
    1543        3143 :     case 3:
    1544        3143 :       A = G3(s2,s3,prec);
    1545        3143 :       if (n == 1) return A;
    1546        2793 :       return divrr(Pi2n(1, prec), mulrr(s3, A));
    1547        1736 :     case 4:
    1548        1736 :       A = G4(s2,prec);
    1549        1736 :       if (n == 1) return A;
    1550        1162 :       return divrr(mulrr(pi, s2), A);
    1551         112 :     case 6:
    1552         112 :       A = sqrr(G3(s2,s3,prec));
    1553         112 :       A = mulrr(A, sqrtr_abs(divsr(3, pi)));
    1554         112 :       A = divrr(A, cbrtu(2, prec));
    1555         112 :       if (n == 1) return A;
    1556          77 :       return divrr(Pi2n(1, prec), A);
    1557          42 :     case 8:
    1558          42 :       A = ellKk(1, s2,s3, prec);
    1559          42 :       B = ellKk(2, s2,s3, prec);
    1560          42 :       A = shiftr(sqrtr_abs(divrr(mulrr(addsr(1, s2), A), sqrtr_abs(pi))), 1);
    1561          42 :       B = shiftr(mulrr(sqrtr_abs(gmul(subrs(s2, 1), mulrr(s2, pi))), B), 3);
    1562          42 :       switch (n)
    1563             :       {
    1564             :         GEN t;
    1565          21 :         case 1: return sqrtr_abs(mulrr(A, B));
    1566           7 :         case 3: return sqrtr_abs(divrr(B, A));
    1567           7 :         case 5: A = sqrtr_abs(divrr(B, A));
    1568           7 :           t = sqrtr_abs(shiftr(addsr(1, shiftr(s2, -1)), -1)); /*sin(3Pi/8)*/
    1569           7 :           return divrr(pi, mulrr(t, A));
    1570           7 :         default: A = sqrtr_abs(mulrr(A, B));
    1571           7 :           t = sqrtr_abs(shiftr(subsr(1, shiftr(s2, -1)), -1)); /*sin(Pi/8)*/
    1572           7 :           return divrr(pi, mulrr(t, A));
    1573             :       }
    1574          42 :     case 12:
    1575          42 :       A = G3(s2,s3,prec);
    1576          42 :       B = G4(s2,prec);
    1577             :       switch (n)
    1578             :       {
    1579             :         GEN t2;
    1580          28 :         case 1: case 11:
    1581          28 :           t2 = shiftr(mulur(27, powrs(divrr(addsr(1,s3), pi), 4)), -2);
    1582          28 :           t2 = mulrr(sqrtnr_abs(t2, 8), mulrr(A, B));
    1583          28 :           if (n == 1) return t2;
    1584           7 :           return divrr(pi, mulrr(sin12(s3), t2));
    1585          14 :         case 5: case 7:
    1586          14 :           t2 = shiftr(divrs(powrs(mulrr(subrs(s3,1), pi), 4), 3), 2);
    1587          14 :           t2 = mulrr(sqrtnr_abs(t2, 8), divrr(B, A));
    1588          14 :           if (n == 5) return t2;
    1589           7 :           return divrr(pi, mulrr(cos12(s3), t2));
    1590             :       }
    1591             :     default: /* n = 24 */
    1592          98 :       if (n > 12)
    1593             :       {
    1594             :         GEN t;
    1595          28 :         n = 24 - n;
    1596          28 :         A = Gn24(n, s2,s3, prec);
    1597             :         switch(n)
    1598             :         { /* t = sin(n*Pi/24) */
    1599           7 :           case 1: t = cos12(s3); t = sinx2(t); break;
    1600           7 :           case 5: t = sin12(s3); t = sinx2(t); break;
    1601           7 :           case 7: t = sin12(s3); togglesign(t); t = sinx2(t); break;
    1602           7 :           default:t = cos12(s3); togglesign(t); t = sinx2(t); break; /* n=11 */
    1603             :         }
    1604          28 :         return divrr(pi, mulrr(A, t));
    1605             :       }
    1606          70 :       return Gn24(n, s2,s3, prec);
    1607             :   }
    1608             : }
    1609             : 
    1610             : /* (a,b) = 1. If 0 < x < b, m >= 0
    1611             : gamma(x/b + m) = gamma(x/b) * mulu_interval_step(x, x+(m-1)*b, b) / b^m
    1612             : gamma(x/b - m) = gamma(x/b) / mulu_interval_step(b-x, b*m-x, b) * (-b)^m */
    1613             : static GEN
    1614       41510 : gammafrac24(GEN a, GEN b, long prec)
    1615             : {
    1616             :   pari_sp av;
    1617             :   long A, B, m, x, bit;
    1618             :   GEN z0, z, t;
    1619       41510 :   if (!(A = itos_or_0(a)) || !(B = itos_or_0(b)) || B > 24) return NULL;
    1620       40572 :   switch(B)
    1621             :   {
    1622       34496 :     case 2: return gammahs(A-1, prec);
    1623        5187 :     case 3: case 4: case 6: case 8: case 12: case 24:
    1624        5187 :       m = A / B;
    1625        5187 :       x = A % B; /* = A - m*B */
    1626        5187 :       if (x < 0) { x += B; m--; } /* now 0 < x < B, A/B = x/B + m */
    1627        5187 :       bit = bit_accuracy(prec);
    1628             :       /* Depending on B and prec, we must experimentally replace the 0.5
    1629             :        * by 0.4 to 2.0 for optimal value. Play safe. */
    1630        5187 :       if (labs(m) > 0.5 * bit * sqrt(bit) / log(bit)) return NULL;
    1631        5173 :       z0 = cgetr(prec); av = avma;
    1632        5173 :       prec += EXTRAPRECWORD;
    1633        5173 :       z = gammafrac24_s(x, B, prec);
    1634        5173 :       if (m)
    1635             :       {
    1636        3619 :         if (m > 0)
    1637        3591 :           t = mpdiv(mulu_interval_step(x, (m-1)*B + x, B), rpowuu(B,m,prec));
    1638             :         else
    1639             :         {
    1640          28 :           m = -m;
    1641          28 :           t = mpdiv(rpowuu(B,m,prec), mulu_interval_step(B-x, m*B - x, B));
    1642          28 :           if (odd(m)) togglesign(t);
    1643             :         }
    1644        3619 :         z = mpmul(z,t);
    1645             :       }
    1646        5173 :       affrr(z, z0); set_avma(av); return z0;
    1647             :   }
    1648         889 :   return NULL;
    1649             : }
    1650             : GEN
    1651      325114 : ggamma(GEN x, long prec)
    1652             : {
    1653             :   pari_sp av;
    1654             :   GEN y;
    1655             : 
    1656      325114 :   switch(typ(x))
    1657             :   {
    1658       61411 :     case t_INT:
    1659       61411 :       if (signe(x) <= 0)
    1660           0 :         pari_err_DOMAIN("gamma","argument", "=",
    1661             :                          strtoGENstr("non-positive integer"), x);
    1662       61411 :       return mpfactr(itos(x) - 1, prec);
    1663             : 
    1664      211373 :     case t_REAL: case t_COMPLEX:
    1665      211373 :       return cxgamma(x, 0, prec);
    1666             : 
    1667       33257 :     case t_FRAC:
    1668             :     {
    1669       33257 :       GEN a = gel(x,1), b = gel(x,2), c = gammafrac24(a, b, prec);
    1670       33257 :       if (c) return c;
    1671        1036 :       av = avma; c = subii(a,b);
    1672        1036 :       if (expi(c) - expi(b) < -50)
    1673             :       { /* x = 1 + c/b is close to 1 */
    1674         161 :         x = mkfrac(c,b);
    1675         161 :         if (lgefint(b) >= prec) x = fractor(x,prec);
    1676         161 :         y = mpexp(lngamma1(x, prec));
    1677             :       }
    1678         875 :       else if (signe(a) < 0 || cmpii(shifti(a,1), b) < 0)
    1679         315 :       { /* gamma will use functional equation x -> z = 1-x = -c/b >= 1/2.
    1680             :          * Gamma(x) = Pi / (sin(Pi z) * Gamma(z)) */
    1681         315 :         GEN z = mkfrac(negi(c), b), q = ground(z), r = gsub(z,q);
    1682         315 :         GEN pi = mppi(prec); /* |r| <= 1/2 */
    1683         315 :         z = fractor(z, prec+EXTRAPRECWORD);
    1684         315 :         y = divrr(pi, mulrr(mpsin(gmul(pi, r)), cxgamma(z, 0, prec)));
    1685         315 :         if (mpodd(q)) togglesign(y);
    1686             :       }
    1687             :       else
    1688             :       {
    1689         560 :         x = fractor(x, prec);
    1690         560 :         y = cxgamma(x, 0, prec);
    1691             :       }
    1692        1036 :       return gerepileupto(av, y);
    1693             :     }
    1694             : 
    1695          70 :     case t_PADIC: return Qp_gamma(x);
    1696       19003 :     default:
    1697       19003 :       av = avma; if (!(y = toser_i(x))) break;
    1698       19005 :       return gerepileupto(av, sergamma(y, prec));
    1699             :   }
    1700           0 :   return trans_eval("gamma",ggamma,x,prec);
    1701             : }
    1702             : 
    1703             : static GEN
    1704         503 : mpfactr_basecase(long n, long prec)
    1705             : {
    1706         503 :   GEN v = cgetg(expu(n) + 2, t_VEC);
    1707         503 :   long k, prec2 = prec + EXTRAPRECWORD;
    1708             :   GEN a;
    1709         503 :   for (k = 1;; k++)
    1710        4213 :   {
    1711        4716 :     long m = n >> (k-1), l;
    1712        4716 :     if (m <= 2) break;
    1713        4213 :     l = (1 + (n >> k)) | 1;
    1714             :     /* product of odd numbers in ]n / 2^k, 2 / 2^(k-1)] */
    1715        4213 :     a = mulu_interval_step_prec(l, m, 2, prec2);
    1716        4213 :     gel(v,k) = k == 1? a: gpowgs(a, k);
    1717             :   }
    1718        4213 :   a = gel(v,--k); while (--k) a = mpmul(a, gel(v,k));
    1719         503 :   if (typ(a) == t_INT) a = itor(a, prec); else a = gprec_wtrunc(a, prec);
    1720         503 :   shiftr_inplace(a, factorial_lval(n, 2));
    1721         503 :   return a;
    1722             : }
    1723             : /* Theory says n > C * b^1.5 / log(b). Timings:
    1724             :  * b = [64, 128, 192, 256, 512, 1024, 2048, 4096, 8192, 16384]
    1725             :  * n = [1930, 2650, 3300, 4270, 9000, 23000, 75000, 210000, 750000, 2400000] */
    1726             : static long
    1727         622 : mpfactr_n(long prec)
    1728             : {
    1729         622 :   long b = bit_accuracy(prec);
    1730         622 :   if (b <=  64) return 1930;
    1731          76 :   if (b <= 128) return 2650;
    1732          55 :   if (b <= 192) return 3300;
    1733          55 :   return b * sqrt(b);
    1734             : }
    1735             : static GEN
    1736        7882 : mpfactr_small(long n, long prec)
    1737             : {
    1738        7882 :   GEN f = cgetr(prec);
    1739        7882 :   pari_sp av = avma;
    1740        7882 :   if (n < 410)
    1741        7882 :     affir(mpfact(n), f);
    1742             :   else
    1743           0 :     affrr(mpfactr_basecase(n, prec), f);
    1744        7882 :   set_avma(av); return f;
    1745             : }
    1746             : GEN
    1747      108986 : mpfactr(long n, long prec)
    1748             : {
    1749      108986 :   GEN f = cgetr(prec);
    1750      108986 :   pari_sp av = avma;
    1751             : 
    1752      108986 :   if (n < 410)
    1753      108364 :     affir(mpfact(n), f);
    1754             :   else
    1755             :   {
    1756         622 :     long N = mpfactr_n(prec);
    1757         503 :     GEN z = n <= N? mpfactr_basecase(n, prec)
    1758         622 :                   : cxgamma(utor(n+1, prec), 0, prec);
    1759         622 :     affrr(z, f);
    1760             :   }
    1761      108986 :   set_avma(av); return f;
    1762             : }
    1763             : 
    1764             : /* First a little worse than mpfactr_n because of the extra logarithm.
    1765             :  * Asymptotically same. */
    1766             : static ulong
    1767        7882 : lngamma_n(long prec)
    1768             : {
    1769        7882 :   long b = bit_accuracy(prec);
    1770             :   double N;
    1771        7882 :   if (b <=  64) return 1450UL;
    1772        7882 :   if (b <= 128) return 2010UL;
    1773         301 :   if (b <= 192) return 2870UL;
    1774         301 :   N = b * sqrt(b);
    1775         301 :   if (b <= 256) return N/1.25;
    1776           0 :   if (b <= 512) return N/1.2;
    1777           0 :   if (b <= 2048) return N/1.1;
    1778           0 :   return N;
    1779             : }
    1780             : 
    1781             : GEN
    1782       35252 : glngamma(GEN x, long prec)
    1783             : {
    1784       35252 :   pari_sp av = avma;
    1785             :   GEN y, y0, t;
    1786             : 
    1787       35252 :   switch(typ(x))
    1788             :   {
    1789        7889 :     case t_INT:
    1790             :     {
    1791             :       ulong n;
    1792        7889 :       if (signe(x) <= 0)
    1793           0 :         pari_err_DOMAIN("lngamma","argument", "=",
    1794             :                          strtoGENstr("non-positive integer"), x);
    1795        7889 :       n = itou_or_0(x);
    1796        7889 :       if (!n || n > lngamma_n(prec)) return cxgamma(x, 1, prec);
    1797        7882 :       return gerepileuptoleaf(av, logr_abs( mpfactr_small(n-1, prec) ));
    1798             :     }
    1799        8253 :     case t_FRAC:
    1800             :     {
    1801        8253 :       GEN a = gel(x,1), b = gel(x,2), c = gammafrac24(a, b, prec);
    1802             :       long e;
    1803        8253 :       if (c) return glog(c, prec);
    1804         805 :       c = subii(a,b); e = expi(b) - expi(c);
    1805         805 :       if (e > 50)
    1806             :       {
    1807           7 :         x = mkfrac(c,b);
    1808           7 :         if (lgefint(b) >= prec) x = fractor(x,prec + nbits2nlong(e));
    1809           7 :         y = lngamma1(x, prec);
    1810             :       }
    1811         798 :       else if (signe(a) < 0 || cmpii(shifti(a,1), b) < 0)
    1812          21 :       { /* gamma will use functional equation x -> z = 1-x = -c/b >= 1/2.
    1813             :          * lngamma(x) = log |Pi / (sin(Pi z) * Gamma(z))| + I*Pi * floor(x) */
    1814          21 :         GEN z = mkfrac(negi(c), b), q = ground(z), r = gsub(z,q);
    1815          21 :         GEN pi = mppi(prec); /* |r| <= 1/2 */
    1816          21 :         z = fractor(z, prec+EXTRAPRECWORD);
    1817          21 :         y = subrr(logr_abs(divrr(pi, mpsin(gmul(pi, r)))), cxgamma(z, 1, prec));
    1818          21 :         if (signe(a) < 0) y = gadd(y, mkcomplex(gen_0, mulri(pi, gfloor(x))));
    1819             :       }
    1820             :       else
    1821             :       {
    1822         777 :         x = fractor(x, e > 1? prec+EXTRAPRECWORD: prec);
    1823         777 :         y = cxgamma(x, 1, prec);
    1824             :       }
    1825         805 :       return gerepileupto(av, y);
    1826             :     }
    1827             : 
    1828       18795 :     case t_REAL: case t_COMPLEX:
    1829       18795 :       return cxgamma(x, 1, prec);
    1830             : 
    1831         308 :     default:
    1832         308 :       if (!(y = toser_i(x))) break;
    1833         308 :       if (lg(y) == 2) pari_err_DOMAIN("lngamma", "argument", "=", gen_0,y);
    1834         301 :       t = serlngamma0(y,prec);
    1835         287 :       y0 = simplify_shallow(gel(y,2));
    1836             :       /* no constant term if y0 = 1 or 2 */
    1837         287 :       if (!isint(y0,&y0) || signe(y0) <= 0 || abscmpiu(y0,2) > 2)
    1838           7 :         t = gadd(t, glngamma(y0,prec));
    1839         287 :       return gerepileupto(av, t);
    1840             : 
    1841           7 :     case t_PADIC: return gerepileupto(av, Qp_log(Qp_gamma(x)));
    1842             :   }
    1843           0 :   return trans_eval("lngamma",glngamma,x,prec);
    1844             : }
    1845             : /********************************************************************/
    1846             : /**                                                                **/
    1847             : /**                  PSI(x) = GAMMA'(x)/GAMMA(x)                   **/
    1848             : /**                                                                **/
    1849             : /********************************************************************/
    1850             : static void
    1851           0 : err_psi(GEN s)
    1852             : {
    1853           0 :   pari_err_DOMAIN("psi","argument", "=",
    1854             :                   strtoGENstr("non-positive integer"), s);
    1855           0 : }
    1856             : static GEN
    1857        4109 : cxpsi(GEN s0, long prec)
    1858             : {
    1859             :   pari_sp av, av2;
    1860             :   GEN sum, z, a, res, sig, tau, s, unr, s2, sq;
    1861             :   long lim, nn, k;
    1862        4109 :   const long la = 3;
    1863        4109 :   int funeq = 0;
    1864             :   pari_timer T;
    1865             : 
    1866        4109 :   if (DEBUGLEVEL>2) timer_start(&T);
    1867        4109 :   s = trans_fix_arg(&prec,&s0,&sig,&tau,&av,&res);
    1868        4109 :   if (signe(sig) <= 0) { funeq = 1; s = gsub(gen_1, s); sig = real_i(s); }
    1869        4109 :   if (typ(s0) == t_INT && signe(s0) <= 0) err_psi(s0);
    1870             : 
    1871        4109 :   if (expo(sig) > 300 || (typ(s) == t_COMPLEX && gexpo(gel(s,2)) > 300))
    1872           7 :   { /* |s| is HUGE. Play safe */
    1873           7 :     GEN L, S = gprec_w(s,LOWDEFAULTPREC), rS = real_i(S), iS = imag_i(S);
    1874             :     double l;
    1875             : 
    1876           7 :     l = rtodbl( gnorm(glog(S, 3)) );
    1877           7 :     l = log(l) / 2.;
    1878           7 :     lim = 2 + (long)ceil((prec2nbits_mul(prec, M_LN2) - l) / (2*(1+log((double)la))));
    1879           7 :     if (lim < 2) lim = 2;
    1880             : 
    1881           7 :     l = (2*lim-1)*la / (2.*M_PI);
    1882           7 :     L = gsub(dbltor(l*l), gsqr(iS));
    1883           7 :     if (signe(L) < 0) L = gen_0;
    1884             : 
    1885           7 :     L = gsub(gsqrt(L, 3), rS);
    1886           7 :     if (signe(L) > 0) nn = (long)ceil(rtodbl(L)); else nn = 1;
    1887           7 :     if (DEBUGLEVEL>2) err_printf("lim, nn: [%ld, %ld]\n",lim,nn);
    1888             :   }
    1889             :   else
    1890             :   {
    1891        4102 :     double ssig = rtodbl(sig);
    1892        4102 :     double st = typ(s) == t_REAL? 0.0: rtodbl(imag_i(s));
    1893             :     double l;
    1894             :     {
    1895             :       double rlog, ilog; /* log (s - Euler) */
    1896        4102 :       dcxlog(ssig - 0.57721566, st, &rlog,&ilog);
    1897        4102 :       l = dnorm(rlog,ilog);
    1898             :     }
    1899        4102 :     if (l < 0.000001) l = 0.000001;
    1900        4102 :     l = log(l) / 2.;
    1901        4102 :     lim = 2 + (long)ceil((prec2nbits_mul(prec, M_LN2) - l) / (2*(1+log((double)la))));
    1902        4102 :     if (lim < 2) lim = 2;
    1903             : 
    1904        4102 :     l = (2*lim-1)*la / (2.*M_PI);
    1905        4102 :     l = l*l - st*st;
    1906        4102 :     if (l < 0.) l = 0.;
    1907        4102 :     nn = (long)ceil( sqrt(l) - ssig );
    1908        4102 :     if (nn < 1) nn = 1;
    1909        4102 :     if (DEBUGLEVEL>2) err_printf("lim, nn: [%ld, %ld]\n",lim,nn);
    1910             :   }
    1911        4109 :   incrprec(prec); unr = real_1(prec); /* one extra word of precision */
    1912        4109 :   s2 = gmul2n(s, 1); sq = gsqr(s);
    1913        4109 :   a = gdiv(unr, gaddgs(s, nn)); /* 1 / (s+n) */
    1914        4109 :   av2 = avma; sum = gmul2n(a, -1);
    1915       97197 :   for (k = 0; k < nn - 1; k += 2)
    1916             :   {
    1917       93088 :     GEN tmp = gaddsg(k*(k + 1), gadd(gmulsg(2*k + 1, s), sq));
    1918       93088 :     sum = gadd(sum, gdiv(gaddsg(2*k + 1, s2), tmp));
    1919       93088 :     if ((k & 1023) == 0) sum = gerepileupto(av2, sum);
    1920             :   }
    1921        4109 :   if (odd(nn)) sum = gadd(sum, gdiv(unr, gaddsg(nn - 1, s)));
    1922        4109 :   z = gsub(glog(gaddgs(s, nn), prec), sum);
    1923        4109 :   if (DEBUGLEVEL>2) timer_printf(&T,"sum from 0 to N - 1");
    1924        4109 :   constbern(lim);
    1925        4109 :   z = gsub(z, psi_sum(gsqr(a), lim));
    1926        4109 :   if (DEBUGLEVEL>2) timer_printf(&T,"Bernoulli sum");
    1927        4109 :   if (funeq)
    1928             :   {
    1929        3990 :     GEN pi = mppi(prec);
    1930        3990 :     z = gadd(z, gmul(pi, gcotan(gmul(pi,s), prec)));
    1931             :   }
    1932        4109 :   set_avma(av); return affc_fixlg(z, res);
    1933             : }
    1934             : 
    1935             : /* n >= 0; return psi(1+x) + O(x^n), x = pol_x(v) */
    1936             : GEN
    1937        9562 : psi1series(long n, long v, long prec)
    1938             : {
    1939        9562 :   long i, l = n+3;
    1940        9562 :   GEN s = cgetg(l, t_SER), z = constzeta(n + 1, prec);
    1941             : 
    1942        9562 :   s[1] = evalsigne(1)|evalvalp(0)|evalvarn(v);
    1943       49854 :   for (i = 1; i <= n+1; i++)
    1944             :   {
    1945       40292 :     GEN c = gel(z,i); /* zeta(i) */
    1946       40292 :     gel(s,i+1) = odd(i)? negr(c): c;
    1947             :   }
    1948        9562 :   return s;
    1949             : }
    1950             : /* T an RgX, return T(X + z0) + O(X^L) */
    1951             : static GEN
    1952     1337302 : tr(GEN T, GEN z0, long L)
    1953             : {
    1954     1337302 :   GEN s = RgX_to_ser(RgX_translate(T, z0), L+3);
    1955     1337302 :   setvarn(s, 0); return s;
    1956             : }
    1957             : /* z0 a complex number with Re(z0) > 1/2; return psi(z0+x) + O(x^L)
    1958             :  * using Luke's rational approximation for psi(x) */
    1959             : static GEN
    1960        5992 : serpsiz0(GEN z0, long L, long v, long prec)
    1961             : {
    1962             :   pari_sp av;
    1963             :   GEN A,A1,A2, B,B1,B2, Q;
    1964             :   long n;
    1965        5992 :   n = gprecision(z0); if (n) prec = n;
    1966        5992 :   z0 = gtofp(z0, prec + EXTRAPRECWORD);
    1967             :   /* Start from n = 3; in Luke's notation, A2 := A_{n-2}, A1 := A_{n-1},
    1968             :    * A := A_n. Same for B */
    1969        5992 :   av = avma;
    1970        5992 :   A2= gdivgs(mkpoln(2, gen_1, utoipos(6)), 2);
    1971        5992 :   B2 = scalarpol_shallow(utoipos(4), 0);
    1972        5992 :   A1= gdivgs(mkpoln(3, gen_1, utoipos(82), utoipos(96)), 6);
    1973        5992 :   B1 = mkpoln(2, utoipos(8), utoipos(28));
    1974        5992 :   A = gdivgs(mkpoln(4, gen_1, utoipos(387), utoipos(2906), utoipos(1920)), 12);
    1975        5992 :   B = mkpoln(3, utoipos(14), utoipos(204), utoipos(310));
    1976        5992 :   A2= tr(A2,z0, L);
    1977        5992 :   B2= tr(B2,z0, L);
    1978        5992 :   A1= tr(A1,z0, L);
    1979        5992 :   B1= tr(B1,z0, L);
    1980        5992 :   A = tr(A, z0, L);
    1981        5992 :   B = tr(B, z0, L); Q = gdiv(A, B);
    1982             :   /* work with z0+x as a variable */
    1983        5992 :   for (n = 4;; n++)
    1984      425794 :   {
    1985      431786 :     GEN Q0 = Q, a, b, r, c3,c2,c1,c0 = muluu(2*n-3, n+1);
    1986      431786 :     GEN u = subiu(muluu(n, 7*n-9), 6);
    1987      431786 :     GEN t = addiu(muluu(n, 7*n-19), 4);
    1988             :     /* c1=(2*n-1)*(3*(n-1)*z+7*n^2-9*n-6);
    1989             :      * c2=(2*n-3)*(z-n-1)*(-3*(n-1)*z+7*n^2-19*n+4);
    1990             :      * c3=(2*n-1)*(n-3)*(z-n)*(z-(n+1))*(z+(n-4)); */
    1991      431786 :     c1 = deg1pol_shallow(muluu(3*(n-1),2*n-1), muliu(u,2*n-1), 0);
    1992      431786 :     c2 = ZX_mul(deg1pol_shallow(utoipos(2*n-3), negi(muluu(2*n-3,n+1)), 0),
    1993      431786 :                 deg1pol_shallow(utoineg(3*(n-1)), t, 0));
    1994      431786 :     r = mkvec3(utoipos(n), utoipos(n+1), stoi(4-n));
    1995      431786 :     c3 = ZX_Z_mul(roots_to_pol(r,0), muluu(2*n-1,n-3));
    1996      431786 :     c1 = tr(c1, z0, L+3);
    1997      431786 :     c2 = tr(c2, z0, L+3);
    1998      431786 :     c3 = tr(c3, z0, L+3);
    1999             : 
    2000             :     /* A_{n+1}, B_{n+1} */
    2001      431786 :     a = gdiv(gadd(gadd(gmul(c1,A),gmul(c2,A1)),gmul(c3,A2)), c0);
    2002      431786 :     b = gdiv(gadd(gadd(gmul(c1,B),gmul(c2,B1)),gmul(c3,B2)), c0);
    2003      431786 :     Q = gdiv(a,b);
    2004      431786 :     if (gexpo(gsub(Q,Q0)) < -prec2nbits(prec)) break;
    2005      425794 :     A2 = A1; A1 = A; A = a;
    2006      425794 :     B2 = B1; B1 = B; B = b;
    2007      425794 :     if (gc_needed(av,1))
    2008             :     {
    2009           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"serpsiz0, n = %ld", n);
    2010           0 :       gerepileall(av, 7, &A,&A1,&A2, &B,&B1,&B2, &Q);
    2011             :     }
    2012             :   }
    2013        5992 :   Q = gmul(Q, gmul2n(gsubsg(1, ginv(tr(pol_x(v),z0, L))), 1));
    2014        5992 :   setvarn(Q, v);
    2015        5992 :   return gadd(negeuler(prec), Q);
    2016             : }
    2017             : /* sum (-1)^k*H(m,k)x^k + O(x^L); L > 0;
    2018             :  * H(m,k) = (-1)^{k * \delta_{m > 0}} sum_{1<=i<m} 1/i^(k+1) */
    2019             : static GEN
    2020        1085 : Hseries(long m, long L, long v, long prec)
    2021             : {
    2022        1085 :   long i, k, bit, l = L+3, M = m < 0? 1-m: m;
    2023        1085 :   pari_sp av = avma;
    2024        1085 :   GEN H = cgetg(l, t_SER);
    2025        1085 :   H[1] = evalsigne(1)|evalvarn(v)|evalvalp(0);
    2026        1085 :   prec++;
    2027        1085 :   bit = -prec2nbits(prec);
    2028        5642 :   for(k = 2; k < l; k++) gel(H,k) = gen_1; /* i=1 */
    2029        1113 :   for (i = 2; i < M; i++)
    2030             :   {
    2031          28 :     GEN ik = invr(utor(i, prec));
    2032         203 :     for (k = 2; k < l; k++)
    2033             :     {
    2034         175 :       if (k > 2) { ik = divru(ik, i); if (expo(ik) < bit) break; }
    2035         175 :       gel(H,k) = gadd(gel(H,k), ik);
    2036             :     }
    2037          28 :     if (gc_needed(av,3))
    2038             :     {
    2039           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"Hseries, i = %ld/%ld", i,M);
    2040           0 :       H = gerepilecopy(av, H);
    2041             :     }
    2042             :   }
    2043        1085 :   if (m > 0)
    2044        3192 :     for (k = 3; k < l; k+=2) togglesign_safe(&gel(H,k));
    2045        1085 :   return H;
    2046             : }
    2047             : 
    2048             : static GEN
    2049       15260 : serpsi(GEN y, long prec)
    2050             : {
    2051       15260 :   GEN Q = NULL, z0, Y = y, Y2;
    2052       15260 :   long L = lg(y)-2, v  = varn(y), vy = valp(y);
    2053             : 
    2054       15260 :   if (!L) pari_err_DOMAIN("psi", "argument", "=", gen_0,y);
    2055       15253 :   if (vy < 0) pari_err_DOMAIN("psi", "series valuation", "<", gen_0,y);
    2056       15253 :   if (vy)
    2057          14 :     z0 = gen_0;
    2058             :   else
    2059             :   {
    2060       15239 :     z0 = simplify_shallow(gel(y,2));
    2061       15239 :     (void)isint(z0, &z0);
    2062             :   }
    2063       15253 :   if (typ(z0) == t_INT && !is_bigint(z0))
    2064             :   {
    2065        9261 :     long m = itos(z0);
    2066        9261 :     if (abscmpiu(muluu(prec2nbits(prec),L), labs(m)) > 0)
    2067             :     { /* psi(m+x) = psi(1+x) + sum_{1 <= i < m} 1/(i+x) for m > 0
    2068             :                     psi(1+x) - sum_{0 <= i < -m} 1/(i+x) for m <= 0 */
    2069        9261 :       GEN H = NULL;
    2070        9261 :       if (m <= 0) L--; /* lose series accuracy due to 1/x term */
    2071        9261 :       if (L)
    2072             :       {
    2073        9254 :         Q = psi1series(L, v, prec);
    2074        9254 :         if (m && m != 1) { H = Hseries(m, L, v, prec); Q = gadd(Q, H); }
    2075        9254 :         if (m <= 0) Q = gsub(Q, ginv(pol_x(v)));
    2076             :       }
    2077             :       else
    2078             :       {
    2079           7 :         Q = scalarser(gen_m1, v, 1);
    2080           7 :         setvalp(Q,-1);
    2081             :       }
    2082             :     }
    2083             :   }
    2084       15253 :   if (!Q)
    2085             :   { /* use psi(1-y)=psi(y)+Pi*cotan(Pi*y) ? */
    2086        5992 :     if (gcmp(real_i(z0),ghalf) < 0) { z0 = gsubsg(1,z0); Y = gsubsg(1,y); }
    2087        5992 :     Q = serpsiz0(z0, L, v, prec);
    2088             :   }
    2089       15253 :   Y2 = serchop0(Y); if (signe(Y2)) Q = gsubst(Q, v, Y2);
    2090             :   /* psi(z0 + Y2) = psi(Y) */
    2091       15253 :   if (Y != y)
    2092             :   { /* psi(y) = psi(Y) + Pi cotan(Pi Y) */
    2093          70 :     GEN pi = mppi(prec);
    2094          70 :     if (typ(z0) == t_INT) Y = Y2; /* in this case cotan(Pi*Y2) = cotan(Pi*Y) */
    2095          70 :     Q = gadd(Q, gmul(pi, gcotan(gmul(pi,Y), prec)));
    2096             :   }
    2097       15253 :   return Q;
    2098             : }
    2099             : 
    2100             : static GEN
    2101        4515 : hrec(ulong a, long b)
    2102             : {
    2103             :   long m;
    2104        4515 :   switch(b - a)
    2105             :   {
    2106        4431 :     case 1: return mkfrac(gen_1, utoipos(a));
    2107          84 :     case 2: return mkfrac(utoipos(2*a + 1), muluu(a, a+1));
    2108             :   }
    2109           0 :   m = (a + b) >> 1;
    2110           0 :   return gadd(hrec(a, m), hrec(m, b));
    2111             : }
    2112             : /* exact Harmonic number H_n */
    2113             : static GEN
    2114        4515 : harmonic(ulong n) { return hrec(1, n+1); }
    2115             : static ulong
    2116       21259 : psi_n(ulong b)
    2117             : {
    2118       21259 :   if (b <= 64) return 50;
    2119       21259 :   if (b <= 128) return 85;
    2120       21259 :   if (b <= 192) return 122;
    2121       21217 :   if (b <= 256) return 150;
    2122       12477 :   if (b <= 512) return 320;
    2123           7 :   if (b <= 1024) return 715;
    2124           0 :   return 0.010709 * pow((double)b, 1.631); /* 1.631 ~ log_3(6) */
    2125             : }
    2126             : GEN
    2127       40761 : gpsi(GEN x, long prec)
    2128             : {
    2129             :   pari_sp av;
    2130             :   ulong n;
    2131             :   GEN y;
    2132       40761 :   switch(typ(x))
    2133             :   {
    2134       21266 :     case t_INT:
    2135       21266 :       if (signe(x) <= 0) err_psi(x);
    2136       21266 :       if (lgefint(x) > 3 || (n = itou(x)) > psi_n(prec2nbits(prec))) break;
    2137       21259 :       av = avma; y = mpeuler(prec);
    2138       21259 :       return gerepileuptoleaf(av, n == 1? negr(y): gsub(harmonic(n-1), y));
    2139        4109 :     case t_REAL: case t_COMPLEX: return cxpsi(x,prec);
    2140       15386 :     default:
    2141       15386 :       av = avma; if (!(y = toser_i(x))) break;
    2142       15260 :       return gerepileupto(av, serpsi(y,prec));
    2143             :   }
    2144         133 :   return trans_eval("psi",gpsi,x,prec);
    2145             : }

Generated by: LCOV version 1.13