Code coverage tests

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

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

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

LCOV - code coverage report
Current view: top level - basemath - trans2.c (source / functions) Hit Total Coverage
Test: PARI/GP v2.8.0 lcov report (development 16912-212c0f0) Lines: 884 916 96.5 %
Date: 2014-10-20 Functions: 57 57 100.0 %
Legend: Lines: hit not hit | Branches: + taken - not taken # not executed Branches: 576 709 81.2 %

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

Generated by: LCOV version 1.9