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

Generated by: LCOV version 1.14