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

Generated by: LCOV version 1.14