Code coverage tests

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

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

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

LCOV - code coverage report
Current view: top level - basemath - trans3.c (source / functions) Hit Total Coverage
Test: PARI/GP v2.18.1 lcov report (development 30316-0578a48613) Lines: 1209 1279 94.5 %
Date: 2025-05-31 09:20:01 Functions: 84 85 98.8 %
Legend: Lines: hit not hit

          Line data    Source code
       1             : /* Copyright (C) 2000  The PARI group.
       2             : 
       3             : This file is part of the PARI/GP package.
       4             : 
       5             : PARI/GP is free software; you can redistribute it and/or modify it under the
       6             : terms of the GNU General Public License as published by the Free Software
       7             : Foundation; either version 2 of the License, or (at your option) any later
       8             : version. It is distributed in the hope that it will be useful, but WITHOUT
       9             : ANY WARRANTY WHATSOEVER.
      10             : 
      11             : Check the License for details. You should have received a copy of it, along
      12             : with the package; see the file 'COPYING'. If not, write to the Free Software
      13             : Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. */
      14             : 
      15             : /********************************************************************/
      16             : /**                                                                **/
      17             : /**                   TRANSCENDENTAL FONCTIONS                     **/
      18             : /**                          (part 3)                              **/
      19             : /**                                                                **/
      20             : /********************************************************************/
      21             : #include "pari.h"
      22             : #include "paripriv.h"
      23             : 
      24             : #define DEBUGLEVEL DEBUGLEVEL_trans
      25             : 
      26             : #define HALF_E 1.3591409 /* exp(1) / 2 */
      27             : 
      28             : /***********************************************************************/
      29             : /**                                                                   **/
      30             : /**                       BESSEL FUNCTIONS                            **/
      31             : /**                                                                   **/
      32             : /***********************************************************************/
      33             : 
      34             : static GEN
      35      900443 : _abs(GEN x)
      36      900443 : { return gabs(gtofp(x,LOWDEFAULTPREC), LOWDEFAULTPREC); }
      37             : /* can we use asymptotic expansion ? */
      38             : static int
      39      408964 : bessel_asymp(GEN n, GEN z, long bit)
      40             : {
      41             :   GEN Z, N;
      42      408964 :   long t = typ(n);
      43      408964 :   if (!is_real_t(t) && t != t_COMPLEX) return 0;
      44      408782 :   Z = _abs(z); N = gaddgs(_abs(n), 1);
      45      408782 :   return gcmpgs(gdiv(Z, gsqr(N)), (bit+10)/2) >= 0; }
      46             : 
      47             : /* Region I: 0 < Arg z <= Pi, II: -Pi < Arg z <= 0 */
      48             : static int
      49         518 : regI(GEN z)
      50             : {
      51         518 :   long s = gsigne(imag_i(z));
      52         518 :   return (s > 0 || (s == 0 && gsigne(real_i(z)) < 0)) ? 1 : 2;
      53             : }
      54             : /* Region 1: Re(z) >= 0, 2: Re(z) < 0, Im(z) >= 0, 3: Re(z) < 0, Im(z) < 0 */
      55             : static int
      56       81899 : regJ(GEN z)
      57             : {
      58       81899 :   if (gsigne(real_i(z)) >= 0) return 1;
      59         336 :   return gsigne(imag_i(z)) >= 0 ? 2 : 3;
      60             : }
      61             : 
      62             : /* Hankel's expansions:
      63             :  * a_k(n) = \prod_{0 <= j < k} (4n^2 - (2j+1)^2)
      64             :  * C(k)[n,z] = a_k(n) / (k! (8 z)^k)
      65             :  * A(z)  = exp(-z) sum_{k >= 0} C(k)
      66             :  * A(-z) = exp(z) sum_{k >= 0} (-1)^k C(k)
      67             :  * J_n(z) ~ [1] (A(z/i) / r + A(-z/i) r) / sqrt(2Pi z)
      68             :  *          [2] (A(z/i) r^3 + A(-z/i) r) / sqrt(2Pi z)
      69             :  *          [3] (A(z/i) / r + A(-z/i) / r^3) / sqrt(2Pi z)
      70             :  * Y_n(z) ~ [1] i(A(z/i) / r + A(-z/i) r) / sqrt(2Pi z)
      71             :  *          [2] i(A(z/i) (r^3-2/r) + A(-z/i) r) / sqrt(2Pi z)
      72             :  *          [3] i(-A(z/i)/r + A(-z/i)(2r-1/r^3)) / sqrt(2Pi z)
      73             :  * K_n(z) ~ A(z) Pi / sqrt(2 Pi z)
      74             :  * I_n(z) ~ [I] (A(-z) + r^2 A(z)) / sqrt(2 Pi z)
      75             :  *          [II](A(-z) + r^(-2) A(z)) / sqrt(2 Pi z) */
      76             : 
      77             : /* set [A(z), A(-z), exp((2*nu+1)*I*Pi/4)] */
      78             : static void
      79       82879 : hankel_ABr(GEN *pA, GEN *pB, GEN *pr, GEN n, GEN z, long bit)
      80             : {
      81       82879 :   GEN E, P, C, Q = gen_0, zi = ginv(gmul2n(z, 3));
      82       82879 :   GEN K = gaddgs(_abs(n), 1), n2 = gmul2n(gsqr(n),2);
      83       82879 :   long prec = nbits2prec(bit), B = bit + 4, m;
      84             : 
      85       82879 :   P = C = real_1_bit(bit);
      86       82879 :   for (m = 1;; m += 2)
      87             :   {
      88     5475740 :     C = gmul(C, gdivgu(gmul(gsub(n2, sqru(2*m - 1)), zi), m));
      89     5475740 :     Q = gadd(Q, C);
      90     5475740 :     C = gmul(C, gdivgu(gmul(gsub(n2, sqru(2*m + 1)), zi), m + 1));
      91     5475740 :     P = gadd(P, C);
      92     5475740 :     if (gexpo(C) < -B && gcmpgs(K, m) <= 0) break;
      93             :   }
      94       82879 :   E = gexp(z, prec);
      95       82879 :   *pA = gdiv(gadd(P, Q), E);
      96       82879 :   *pB = gmul(gsub(P, Q), E);
      97       82879 :   *pr = gexp(mulcxI(gmul(gaddgs(gmul2n(n,1), 1), Pi2n(-2, prec))), prec);
      98       82879 : }
      99             : 
     100             : /* sqrt(2*Pi*z) */
     101             : static GEN
     102       82879 : sqz(GEN z, long bit)
     103             : {
     104       82879 :   long prec = nbits2prec(bit);
     105       82879 :   return gsqrt(gmul(Pi2n(1, prec), z), prec);
     106             : }
     107             : 
     108             : static GEN
     109         462 : besskasymp(GEN nu, GEN z, long bit)
     110             : {
     111             :   GEN A, B, r;
     112         462 :   long prec = nbits2prec(bit);
     113         462 :   hankel_ABr(&A,&B,&r, nu, z, bit);
     114         462 :   return gdiv(gmul(A, mppi(prec)), sqz(z, bit));
     115             : }
     116             : 
     117             : static GEN
     118         518 : bessiasymp(GEN nu, GEN z, long bit)
     119             : {
     120             :   GEN A, B, r, R, r2;
     121         518 :   hankel_ABr(&A,&B,&r, nu, z, bit);
     122         518 :   r2 = gsqr(r);
     123         518 :   R = regI(z) == 1 ? gmul(A, r2) : gdiv(A, r2);
     124         518 :   return gdiv(gadd(B, R), sqz(z, bit));
     125             : }
     126             : 
     127             : static GEN
     128       81433 : bessjasymp(GEN nu, GEN z, long bit)
     129             : {
     130             :   GEN A, B, r, R;
     131       81433 :   long reg = regJ(z);
     132       81433 :   hankel_ABr(&A,&B,&r, nu, mulcxmI(z), bit);
     133       81433 :   if (reg == 1) R = gadd(gdiv(A, r), gmul(B, r));
     134         168 :   else if (reg == 2) R = gadd(gmul(A, gpowgs(r, 3)), gmul(B, r));
     135          56 :   else R = gadd(gdiv(A, r), gdiv(B, gpowgs(r, 3)));
     136       81433 :   return gdiv(R, sqz(z, bit));
     137             : }
     138             : 
     139             : static GEN
     140         466 : bessyasymp(GEN nu, GEN z, long bit)
     141             : {
     142             :   GEN A, B, r, R;
     143         466 :   long reg = regJ(z);
     144         466 :   hankel_ABr(&A,&B,&r, nu, mulcxmI(z), bit);
     145         466 :   if (reg == 1) R = gsub(gmul(B, r), gdiv(A, r));
     146         168 :   else if (reg == 2)
     147         112 :     R = gadd(gmul(A, gsub(gpowgs(r, 3), gmul2n(ginv(r), 1))), gmul(B, r));
     148             :   else
     149          56 :     R = gsub(gmul(B, gsub(gmul2n(r, 1), ginv(gpowgs(r, 3)))), gdiv(A, r));
     150         466 :   return gdiv(mulcxI(R), sqz(z, bit));
     151             : }
     152             : 
     153             : /* n! sum_{0 <= k <= m} x^k / (k!*(k+n)!) */
     154             : static GEN
     155      314931 : _jbessel(GEN n, GEN x, long m)
     156             : {
     157      314931 :   pari_sp av = avma;
     158      314931 :   GEN s = gen_1;
     159             :   long k;
     160             : 
     161   109563028 :   for (k = m; k >= 1; k--)
     162             :   {
     163   109248097 :     s = gaddsg(1, gdiv(gmul(x,s), gmulgu(gaddgs(n, k), k)));
     164   109248097 :     if (gc_needed(av,1))
     165             :     {
     166           0 :       if (DEBUGMEM>1) pari_warn(warnmem,"besselj");
     167           0 :       s = gc_upto(av, s);
     168             :     }
     169             :   }
     170      314931 :   return s;
     171             : }
     172             : 
     173             : /* max(2, L * approximate solution to x log x = B) */
     174             : static long
     175      325504 : bessel_get_lim(double B, double L)
     176      325504 : { return maxss(2, L * exp(dbllambertW0(B))); }
     177             : 
     178          42 : static GEN vjbesselh(void* E, GEN z, long prec){return jbesselh((GEN)E,z,prec);}
     179         126 : static GEN vjbessel(void* E, GEN z, long prec) {return jbessel((GEN)E,z,prec);}
     180          42 : static GEN vibessel(void* E, GEN z, long prec) {return ibessel((GEN)E,z,prec);}
     181         126 : static GEN vnbessel(void* E, GEN z, long prec) {return nbessel((GEN)E,z,prec);}
     182          42 : static GEN vkbessel(void* E, GEN z, long prec) {return kbessel((GEN)E,z,prec);}
     183             : 
     184             : /* if J != 0 BesselJ, else BesselI. */
     185             : static GEN
     186      397197 : jbesselintern(GEN n, GEN z, long J, long prec)
     187             : {
     188      397197 :   const char *f = J? "besselj": "besseli";
     189             :   long i, ki;
     190      397197 :   pari_sp av = avma;
     191             :   GEN y;
     192             : 
     193      397197 :   switch(typ(z))
     194             :   {
     195      396791 :     case t_INT: case t_FRAC: case t_REAL: case t_COMPLEX:
     196             :     {
     197      396791 :       int flz0 = gequal0(z);
     198             :       long lim, k, precnew, bit;
     199             :       GEN p1, p2;
     200             :       double az, L;
     201             : 
     202      396791 :       i = precision(z); if (i) prec = i;
     203      396791 :       if (flz0 && gequal0(n)) return real_1(prec);
     204      396791 :       bit = prec2nbits(prec);
     205      396791 :       if (bessel_asymp(n, z, bit))
     206             :       {
     207       81951 :         GEN R = J? bessjasymp(n, z, bit): bessiasymp(n, z, bit);
     208       81951 :         if (typ(R) == t_COMPLEX && isexactzero(imag_i(n))
     209       81265 :                                 && gsigne(real_i(z)) > 0
     210       81111 :                                 && isexactzero(imag_i(z))) R = gcopy(gel(R,1));
     211       81951 :         return gc_upto(av, R);
     212             :       }
     213      314840 :       p2 = gpow(gmul2n(z,-1),n,prec);
     214      314812 :       p2 = gdiv(p2, ggamma(gaddgs(n,1),prec));
     215      314812 :       if (flz0) return gc_upto(av, p2);
     216      314812 :       az = dblmodulus(z); L = HALF_E * az;
     217      314812 :       precnew = prec;
     218      314812 :       if (az >= 1.0) precnew += 1 + nbits2extraprec((long)(az/M_LN2));
     219      314812 :       if (issmall(n,&ki)) {
     220      313881 :         k = labs(ki);
     221      313881 :         n = utoi(k);
     222             :       } else {
     223         847 :         i = precision(n);
     224         847 :         if (i && i < precnew) n = gtofp(n,precnew);
     225             :       }
     226      314728 :       z = gtofp(z,precnew);
     227      314728 :       lim = bessel_get_lim(prec2nbits_mul(prec,M_LN2/2) / L, L);
     228      314728 :       z = gmul2n(gsqr(z),-2); if (J) z = gneg(z);
     229      314728 :       p1 = gprec_wtrunc(_jbessel(n,z,lim), prec);
     230      314728 :       return gc_upto(av, gmul(p2,p1));
     231             :     }
     232             : 
     233          14 :     case t_PADIC: pari_err_IMPL(stack_strcat("p-adic ",f));
     234         392 :     default:
     235             :     {
     236             :       long v, k, m;
     237         392 :       if (!(y = toser_i(z))) break;
     238         238 :       if (issmall(n,&ki)) n = utoi(labs(ki));
     239         210 :       y = gmul2n(gsqr(y),-2); if (J) y = gneg(y);
     240         210 :       v = valser(y);
     241         210 :       if (v < 0) pari_err_DOMAIN(f, "valuation", "<", gen_0, z);
     242         203 :       if (v == 0) pari_err_IMPL(stack_strcat(f, " around a!=0"));
     243         203 :       m = lg(y) - 2;
     244         203 :       k = m - (v >> 1);
     245         203 :       if (k <= 0) { set_avma(av); return scalarser(gen_1, varn(z), v); }
     246         203 :       setlg(y, k+2); return gc_upto(av, _jbessel(n, y, m));
     247             :     }
     248             :   }
     249         154 :   return trans_evalgen(f, (void*)n, J? vjbessel: vibessel, z, prec);
     250             : }
     251             : GEN
     252      384972 : jbessel(GEN n, GEN z, long prec) { return jbesselintern(n,z,1,prec); }
     253             : GEN
     254         896 : ibessel(GEN n, GEN z, long prec) { return jbesselintern(n,z,0,prec); }
     255             : 
     256             : /* k > 0 */
     257             : static GEN
     258         119 : _jbesselh(long k, GEN z, long prec)
     259             : {
     260         119 :   GEN s, c, p0, p1, zinv = ginv(z);
     261             :   long i;
     262             : 
     263         119 :   gsincos(z,&s,&c,prec);
     264         119 :   p1 = gmul(zinv,s);
     265         119 :   p0 = p1; p1 = gmul(zinv,gsub(p0,c));
     266        1134 :   for (i = 2; i <= k; i++)
     267             :   {
     268        1015 :     GEN p2 = gsub(gmul(gmulsg(2*i-1,zinv), p1), p0);
     269        1015 :     p0 = p1; p1 = p2;
     270             :   }
     271         119 :   return p1;
     272             : }
     273             : 
     274             : /* J_{n+1/2}(z) */
     275             : GEN
     276         315 : jbesselh(GEN n, GEN z, long prec)
     277             : {
     278             :   long k, i;
     279             :   pari_sp av;
     280             :   GEN y;
     281             : 
     282         315 :   if (typ(n)!=t_INT) pari_err_TYPE("jbesselh",n);
     283         203 :   k = itos(n);
     284         203 :   if (k < 0) return jbessel(gadd(ghalf,n), z, prec);
     285             : 
     286         203 :   switch(typ(z))
     287             :   {
     288         133 :     case t_INT: case t_FRAC: case t_REAL: case t_COMPLEX:
     289             :     {
     290             :       long pr;
     291             :       GEN p1;
     292         133 :       if (gequal0(z))
     293             :       {
     294           7 :         av = avma;
     295           7 :         p1 = gmul(gsqrt(gdiv(z,mppi(prec)),prec),gpowgs(z,k));
     296           7 :         p1 = gdiv(p1, mulu_interval(k+1, 2*k+1)); /* x k! / (2k+1)! */
     297           7 :         return gc_upto(av, gmul2n(p1,2*k));
     298             :       }
     299         126 :       if ( (pr = precision(z)) ) prec = pr;
     300         126 :       if (bessel_asymp(n, z, prec2nbits(prec)))
     301           7 :         return jbessel(gadd(ghalf,n), z, prec);
     302         119 :       y = cgetc(prec); av = avma;
     303         119 :       p1 = gsqrt(gdiv(z, Pi2n(-1,prec)), prec);
     304         119 :       if (!k)
     305          21 :         p1 = gmul(p1, gsinc(z, prec));
     306             :       else
     307             :       {
     308          98 :         long bits = BITS_IN_LONG + 2*k * (log2(k) -  dbllog2(z));
     309          98 :         if (bits > 0)
     310             :         {
     311          98 :           prec += nbits2extraprec(bits);
     312          98 :           if (pr) z = gtofp(z, prec);
     313             :         }
     314          98 :         p1 = gmul(p1, _jbesselh(k,z,prec));
     315             :       }
     316         119 :       set_avma(av); return affc_fixlg(p1, y);
     317             :     }
     318           0 :     case t_PADIC: pari_err_IMPL("p-adic jbesselh function");
     319          70 :     default:
     320             :     {
     321             :       long t, v;
     322          70 :       av = avma; if (!(y = toser_i(z))) break;
     323          35 :       if (gequal0(y)) return gc_upto(av, gpowgs(y,k));
     324          35 :       v = valser(y);
     325          35 :       if (v < 0) pari_err_DOMAIN("besseljh","valuation", "<", gen_0, z);
     326          28 :       t = lg(y)-2;
     327          28 :       if (v) y = sertoser(y, t + (2*k+1)*v);
     328          28 :       if (!k)
     329           7 :         y = gsinc(y,prec);
     330             :       else
     331             :       {
     332          21 :         GEN T, a = _jbesselh(k, y, prec);
     333          21 :         if (v) y = sertoser(y, t + k*v); /* lower precision */
     334          21 :         y = gdiv(a, gpowgs(y, k));
     335          21 :         T = cgetg(k+1, t_VECSMALL);
     336         168 :         for (i = 1; i <= k; i++) T[i] = 2*i+1;
     337          21 :         y = gmul(y, zv_prod_Z(T));
     338             :       }
     339          28 :       return gc_upto(av, y);
     340             :     }
     341             :   }
     342          35 :   return trans_evalgen("besseljh",(void*)n, vjbesselh, z, prec);
     343             : }
     344             : 
     345             : static GEN
     346           0 : kbessel2(GEN nu, GEN x, long prec)
     347             : {
     348           0 :   pari_sp av = avma;
     349           0 :   GEN p1, a, x2 = gshift(x,1);
     350             : 
     351           0 :   a = gtofp(gaddgs(gshift(nu,1), 1), prec);
     352           0 :   p1 = hyperu(gshift(a,-1), a, x2, prec);
     353           0 :   p1 = gmul(gmul(p1, gpow(x2,nu,prec)), sqrtr(mppi(prec)));
     354           0 :   return gc_upto(av, gmul(p1, gexp(gneg(x),prec)));
     355             : }
     356             : 
     357             : /* special case of hyperu */
     358             : static GEN
     359          14 : kbessel1(GEN nu, GEN gx, long prec)
     360             : {
     361             :   GEN x, y, zf, r, u, pi, nu2;
     362          14 :   long bit, k, k2, n2, n, l = (typ(gx)==t_REAL)? realprec(gx): prec;
     363             :   pari_sp av;
     364             : 
     365          14 :   if (typ(nu)==t_COMPLEX) return kbessel2(nu, gx, l);
     366          14 :   y = cgetr(l); av = avma;
     367          14 :   x = gtofp(gx, l);
     368          14 :   nu = gtofp(nu,l); nu2 = sqrr(nu);
     369          14 :   shiftr_inplace(nu2,2); togglesign(nu2); /* nu2 = -4nu^2 */
     370          14 :   n = (long) (prec2nbits_mul(l,M_LN2) + M_PI*fabs(rtodbl(nu))) / 2;
     371          14 :   bit = prec2nbits(l) - 1;
     372          14 :   l += EXTRAPREC64;
     373          14 :   pi = mppi(l); n2 = n<<1; r = gmul2n(x,1);
     374          14 :   if (cmprs(x, n) < 0)
     375             :   {
     376          14 :     pari_sp av2 = avma;
     377          14 :     GEN q, v, c, s = real_1(l), t = real_0(l);
     378        1246 :     for (k = n2, k2 = 2*n2-1; k > 0; k--, k2 -= 2)
     379             :     {
     380        1232 :       GEN ak = divri(addri(nu2, sqru(k2)), mulss(n2<<2, -k));
     381        1232 :       s = addsr(1, mulrr(ak,s));
     382        1232 :       t = addsr(k2,mulrr(ak,t));
     383        1232 :       if (gc_needed(av2,3)) (void)gc_all(av2, 2, &s,&t);
     384             :     }
     385          14 :     shiftr_inplace(t, -1);
     386          14 :     q = utor(n2, l);
     387          14 :     zf = sqrtr(divru(pi,n2));
     388          14 :     u = gprec_wensure(mulrr(zf, s), l);
     389          14 :     v = gprec_wensure(divrs(addrr(mulrr(t,zf),mulrr(u,nu)),-n2), l);
     390             :     for(;;)
     391         301 :     {
     392         315 :       GEN p1, e, f, d = real_1(l);
     393             :       pari_sp av3;
     394         315 :       c = divur(5,q); if (expo(c) >= -1) c = real2n(-1,l);
     395         315 :       p1 = subsr(1, divrr(r,q)); if (cmprr(c,p1)>0) c = p1;
     396         315 :       togglesign(c); av3 = avma;
     397         315 :       e = u;
     398         315 :       f = v;
     399         315 :       for (k = 1;; k++)
     400       35230 :       {
     401       35545 :         GEN w = addrr(gmul2n(mulur(2*k-1,u), -1), mulrr(subrs(q,k),v));
     402       35545 :         w = addrr(w, mulrr(nu, subrr(u,gmul2n(v,1))));
     403       35545 :         u = divru(mulrr(q,v), k);
     404       35545 :         v = divru(w,k);
     405       35545 :         d = mulrr(d,c);
     406       35545 :         e = addrr(e, mulrr(d,u));
     407       35545 :         f = addrr(f, p1 = mulrr(d,v));
     408       35545 :         if (expo(p1) - expo(f) <= 1-prec2nbits(realprec(p1))) break;
     409       35230 :         if (gc_needed(av3,3)) (void)gc_all(av3,5,&u,&v,&d,&e,&f);
     410             :       }
     411         315 :       u = e;
     412         315 :       v = f;
     413         315 :       q = mulrr(q, addrs(c,1));
     414         315 :       if (expo(r) - expo(subrr(q,r)) >= bit) break;
     415         301 :       (void)gc_all(av2, 3, &u,&v,&q);
     416             :     }
     417          14 :     u = mulrr(u, gpow(divru(x,n),nu,prec));
     418             :   }
     419             :   else
     420             :   {
     421           0 :     GEN s, zz = ginv(gmul2n(r,2));
     422           0 :     pari_sp av2 = avma;
     423           0 :     s = real_1(l);
     424           0 :     for (k = n2, k2 = 2*n2-1; k > 0; k--, k2 -= 2)
     425             :     {
     426           0 :       GEN ak = divru(mulrr(addri(nu2, sqru(k2)), zz), k);
     427           0 :       s = subsr(1, mulrr(ak,s));
     428           0 :       if (gc_needed(av2,3)) s = gc_leaf(av2, s);
     429             :     }
     430           0 :     zf = sqrtr(divrr(pi,r));
     431           0 :     u = mulrr(s, zf);
     432             :   }
     433          14 :   affrr(mulrr(u, mpexp(negr(x))), y);
     434          14 :   set_avma(av); return y;
     435             : }
     436             : 
     437             : /*   sum_{k=0}^m x^k (H(k)+H(k+n)) / (k! (k+n)!)
     438             :  * + sum_{k=0}^{n-1} (-x)^(k-n) (n-k-1)!/k! */
     439             : static GEN
     440       10860 : _kbessel(long n, GEN x, long m, long prec)
     441             : {
     442             :   GEN p1, p2, s, H;
     443       10860 :   long k, M = m + n, exact = (M <= prec2nbits(prec));
     444             :   pari_sp av;
     445             : 
     446       10860 :   H = cgetg(M+2,t_VEC); gel(H,1) = gen_0;
     447       10860 :   if (exact)
     448             :   {
     449       10853 :     gel(H,2) = s = gen_1;
     450      479307 :     for (k=2; k<=M; k++) gel(H,k+1) = s = gdivgu(gaddsg(1,gmulsg(k,s)),k);
     451             :   }
     452             :   else
     453             :   {
     454           7 :     gel(H,2) = s = real_1(prec);
     455        2877 :     for (k=2; k<=M; k++) gel(H,k+1) = s = divru(addsr(1,mulur(k,s)),k);
     456             :   }
     457       10860 :   s = gadd(gel(H,m+1), gel(H,M+1)); av = avma;
     458      430240 :   for (k = m; k > 0; k--)
     459             :   {
     460      419380 :     s = gadd(gadd(gel(H,k),gel(H,k+n)), gdiv(gmul(x,s),mulss(k,k+n)));
     461      419380 :     if (gc_needed(av,1))
     462             :     {
     463           0 :       if (DEBUGMEM>1) pari_warn(warnmem,"_kbessel");
     464           0 :       s = gc_upto(av, s);
     465             :     }
     466             :   }
     467       10860 :   p1 = exact? mpfact(n): mpfactr(n,prec);
     468       10860 :   s = gdiv(s,p1);
     469       10860 :   if (n)
     470             :   {
     471        8330 :     x = gneg(ginv(x));
     472        8330 :     p2 = gmulsg(n, gdiv(x,p1));
     473        8330 :     s = gadd(s,p2);
     474       62804 :     for (k=n-1; k>0; k--)
     475             :     {
     476       54474 :       p2 = gmul(p2, gmul(mulss(k,n-k),x));
     477       54474 :       s = gadd(s,p2);
     478             :     }
     479             :   }
     480       10860 :   return s;
     481             : }
     482             : 
     483             : /* N = 1: Bessel N, else Bessel K */
     484             : static GEN
     485       12432 : kbesselintern(GEN n, GEN z, long N, long prec)
     486             : {
     487       12432 :   const char *f = N? "besseln": "besselk";
     488             :   long i, k, ki, lim, precnew, fl2, ex, bit;
     489       12432 :   pari_sp av = avma;
     490             :   GEN p1, p2, y, p3, pp, pm, s, c;
     491             :   double az;
     492             : 
     493       12432 :   switch(typ(z))
     494             :   {
     495       12047 :     case t_INT: case t_FRAC: case t_REAL: case t_COMPLEX:
     496       12047 :       if (gequal0(z)) pari_err_DOMAIN(f, "argument", "=", gen_0, z);
     497       12047 :       i = precision(z); if (i) prec = i;
     498       12047 :       i = precision(n); if (i && prec > i) prec = i;
     499       12047 :       bit = prec2nbits(prec);
     500       12047 :       if (bessel_asymp(n, z, bit))
     501             :       {
     502         928 :         GEN R = N? bessyasymp(n, z, bit): besskasymp(n, z, bit);
     503         928 :         if (typ(R) == t_COMPLEX && isexactzero(imag_i(n))
     504         221 :                                 && gsigne(real_i(z)) > 0
     505          81 :                                 && isexactzero(imag_i(z))) R = gcopy(gel(R,1));
     506         928 :         return gc_upto(av, R);
     507             :       }
     508             :       /* heuristic threshold */
     509       11119 :       if (!N && !gequal0(n) && gexpo(z) > bit/16 + gexpo(n))
     510          14 :         return kbessel1(n,z,prec);
     511       11098 :       az = dblmodulus(z); precnew = prec;
     512       11098 :       if (az >= 1) precnew += 1 + nbits2extraprec((long)((N?az:2*az)/M_LN2));
     513       11098 :       z = gtofp(z, precnew);
     514       11098 :       if (issmall(n,&ki))
     515             :       {
     516       10776 :         GEN z2 = gmul2n(z, -1), Z;
     517       10776 :         double B, L = HALF_E * az;
     518       10776 :         k = labs(ki);
     519       10776 :         B = prec2nbits_mul(prec,M_LN2/2) / L;
     520       10776 :         if (!N) B += 0.367879; /* exp(-1) */
     521       10776 :         lim = bessel_get_lim(B, L);
     522       10776 :         Z = gsqr(z2); if (N) Z = gneg(Z);
     523       10776 :         p1 = gmul(gpowgs(z2,k), _kbessel(k, Z, lim, precnew));
     524       10776 :         p2 = gadd(mpeuler(precnew), glog(z2,precnew));
     525       10776 :         p3 = jbesselintern(stoi(k),z,N,precnew);
     526       10776 :         p2 = gsub(gmul2n(p1,-1),gmul(p2,p3));
     527       10776 :         p2 = gprec_wtrunc(p2, prec);
     528       10776 :         if (N)
     529             :         {
     530        8361 :           p2 = gdiv(p2, Pi2n(-1,prec));
     531        8361 :           if (ki >= 0 || !odd(k)) p2 = gneg(p2);
     532             :         } else
     533        2415 :           if (odd(k)) p2 = gneg(p2);
     534       10776 :         return gc_GEN(av, p2);
     535             :       }
     536             : 
     537         259 :       n = gtofp(n, precnew);
     538         259 :       gsincos(gmul(n,mppi(precnew)), &s,&c,precnew);
     539         259 :       ex = gexpo(s);
     540         259 :       if (ex < 0) precnew += nbits2extraprec(N? -ex: -2*ex);
     541         259 :       if (i && i < precnew) {
     542          84 :         n = gtofp(n,precnew);
     543          84 :         z = gtofp(z,precnew);
     544          84 :         gsincos(gmul(n,mppi(precnew)), &s,&c,precnew);
     545             :       }
     546             : 
     547         259 :       pp = jbesselintern(n,      z,N,precnew);
     548         259 :       pm = jbesselintern(gneg(n),z,N,precnew);
     549         259 :       if (N)
     550         189 :         p1 = gsub(gmul(c,pp),pm);
     551             :       else
     552          70 :         p1 = gmul(gsub(pm,pp), Pi2n(-1,precnew));
     553         259 :       p1 = gdiv(p1, s);
     554         259 :       return gc_GEN(av, gprec_wtrunc(p1,prec));
     555             : 
     556          14 :     case t_PADIC: pari_err_IMPL(stack_strcat("p-adic ",f));
     557         371 :     default:
     558         371 :       if (!(y = toser_i(z))) break;
     559         217 :       if (issmall(n,&ki))
     560             :       {
     561         105 :         long v, mv, k = labs(ki), m = lg(y)-2;
     562         105 :         y = gmul2n(gsqr(y),-2); if (N) y = gneg(y);
     563         105 :         v = valser(y);
     564         105 :         if (v < 0) pari_err_DOMAIN(f, "valuation", "<", gen_0, z);
     565          91 :         if (v == 0) pari_err_IMPL(stack_strcat(f, " around a!=0"));
     566          91 :         mv = m - (v >> 1);
     567          91 :         if (mv <= 0) { set_avma(av); return scalarser(gen_1, varn(z), v); }
     568          84 :         setlg(y, mv+2); return gc_GEN(av, _kbessel(k, y, m, prec));
     569             :       }
     570          98 :       if (!issmall(gmul2n(n,1),&ki))
     571          70 :         pari_err_DOMAIN(f, "2n mod Z", "!=", gen_0, n);
     572          28 :       k = labs(ki); n = gmul2n(stoi(k),-1);
     573          28 :       fl2 = (k&3)==1;
     574          28 :       pm = jbesselintern(gneg(n), y, N, prec);
     575          28 :       if (N) p1 = pm;
     576             :       else
     577             :       {
     578           7 :         pp = jbesselintern(n, y, N, prec);
     579           7 :         p2 = gpowgs(y,-k); if (fl2 == 0) p2 = gneg(p2);
     580           7 :         p3 = gmul2n(diviiexact(mpfact(k + 1),mpfact((k + 1) >> 1)),-(k + 1));
     581           7 :         p3 = gdivgu(gmul2n(gsqr(p3),1),k);
     582           7 :         p2 = gmul(p2,p3);
     583           7 :         p1 = gsub(pp,gmul(p2,pm));
     584             :       }
     585          28 :       return gc_upto(av, fl2? gneg(p1): gcopy(p1));
     586             :   }
     587         154 :   return trans_evalgen(f, (void*)n, N? vnbessel: vkbessel, z, prec);
     588             : }
     589             : 
     590             : GEN
     591        3115 : kbessel(GEN n, GEN z, long prec) { return kbesselintern(n,z,0,prec); }
     592             : GEN
     593        9317 : ybessel(GEN n, GEN z, long prec) { return kbesselintern(n,z,1,prec); }
     594             : /* J + iN */
     595             : GEN
     596         224 : hbessel1(GEN n, GEN z, long prec)
     597             : {
     598         224 :   pari_sp av = avma;
     599         224 :   GEN J = jbessel(n,z,prec);
     600         196 :   GEN Y = ybessel(n,z,prec);
     601         182 :   return gc_upto(av, gadd(J, mulcxI(Y)));
     602             : }
     603             : /* J - iN */
     604             : GEN
     605         224 : hbessel2(GEN n, GEN z, long prec)
     606             : {
     607         224 :   pari_sp av = avma;
     608         224 :   GEN J = jbessel(n,z,prec);
     609         196 :   GEN Y = ybessel(n,z,prec);
     610         182 :   return gc_upto(av, gadd(J, mulcxmI(Y)));
     611             : }
     612             : 
     613             : static GEN
     614        1008 : besselrefine(GEN z, GEN nu, GEN (*B)(GEN,GEN,long), long bit)
     615             : {
     616        1008 :   GEN z0 = gprec_w(z, DEFAULTPREC), nu1 = gaddgs(nu, 1), t;
     617        1008 :   long e, n, c, j, prec = DEFAULTPREC;
     618             : 
     619        1008 :   t = gdiv(B(nu1, z0, prec), B(nu, z0, prec));
     620        1008 :   t = gadd(z0, gdiv(gsub(gsqr(z0), gsqr(nu)), gsub(gdiv(nu, z0), t)));
     621        1008 :   e = gexpo(t) - 2 * gexpo(z0) - 1; if (e < 0) e = 0;
     622        1008 :   n = expu(bit + 32 - e);
     623        1008 :   c = 1 + e + ((bit - e) >> n);
     624        8064 :   for (j = 1; j <= n; j++)
     625             :   {
     626        7056 :     c = 2 * c - e;
     627        7056 :     prec = nbits2prec(c); z = gprec_w(z, prec);
     628        7056 :     t = gdiv(B(nu1, z, prec), B(nu, z, prec));
     629        7056 :     z = gsub(z, ginv(gsub(gdiv(nu, z), t)));
     630             :   }
     631        1008 :   return gprec_w(z, nbits2prec(bit));
     632             : }
     633             : 
     634             : /* solve tan(fi) - fi = y, y >= 0; Temme's method */
     635             : static double
     636         700 : fi(double y)
     637             : {
     638             :   double p, pp, r;
     639         700 :   if (y == 0) return 0;
     640         700 :   if (y > 100000) return M_PI/2;
     641         700 :   if (y < 1)
     642             :   {
     643         455 :     p = pow(3*y, 1.0/3); pp = p * p;
     644         455 :     p = p * (1 + pp * (-210 * pp + (27 - 2*pp)) / 1575);
     645             :   }
     646             :   else
     647             :   {
     648         245 :     p = 1 / (y + M_PI/2); pp = p * p;
     649         245 :     p = M_PI/2 - p*(1 + pp*(2310 + pp*(3003 + pp*(4818 + pp*(8591 + pp*16328)))) / 3465);
     650             :   }
     651         700 :   pp = (y + p) * (y + p); r = (p - atan(p + y)) / pp;
     652         700 :   return p - (1 + pp) * r * (1 + r / (p + y));
     653             : }
     654             : 
     655             : static GEN
     656        1022 : besselzero(GEN nu, long n, GEN (*B)(GEN,GEN,long), long bit)
     657             : {
     658        1022 :   pari_sp av = avma;
     659        1022 :   long prec = nbits2prec(bit);
     660        1022 :   int J = B == jbessel;
     661             :   GEN z;
     662        1022 :   if (n <= 0) pari_err_DOMAIN("besselzero", "n", "<=", gen_0, stoi(n));
     663        1008 :   if (n > LONG_MAX / 4) pari_err_OVERFLOW("besselzero");
     664        1008 :   if (is_real_t(typ(nu)) && gsigne(nu) >= 0)
     665        1008 :   { /* Temme */
     666        1008 :     double x, c, b, a = gtodouble(nu), t = J? 0.25: 0.75;
     667        1008 :     if (n >= 3*a - 8)
     668             :     {
     669         308 :       double aa = a*a, mu = 4*aa, mu2 = mu*mu, p, p0, p1, q1;
     670         308 :       p = 7 * mu - 31; p0 = mu-1;
     671         308 :       if (1 + p == p) /* p large */
     672           0 :         p1 = q1 = 0;
     673             :       else
     674             :       {
     675         308 :         p1 = 4 * (253 * mu2 - 3722 * mu + 17869) / (15 * p);
     676         308 :         q1 = 1.6 * (83 * mu2 - 982 * mu + 3779) / p;
     677             :       }
     678         308 :       b = (n + a/2 - t) * M_PI;
     679         308 :       c = 1 / (64 * b * b);
     680         308 :       x = b - p0 * (1 - p1 * c) / (8 * b * (1 - q1 * c));
     681             :     }
     682             :     else
     683             :     {
     684         700 :       double u, v, w, xx, bb = a >= 3? pow(a, -2./3): 1;
     685         700 :       if (n == 1)
     686         336 :         x = J? -2.33811: -1.17371;
     687             :       else
     688             :       {
     689         364 :         double pp1 = 5./48, qq1 = -5./36, y = 3./8 * M_PI;
     690         364 :         x = 4 * y * (n - t); v = 1 / (x*x);
     691         364 :         x = - pow(x, 2.0/3) * (1 + v * (pp1 + qq1 * v));
     692             :       }
     693         700 :       u = x * bb; v = fi(2.0/3 * pow(-u, 1.5));
     694         700 :       w = 1 / cos(v); xx = 1 - w*w; c = sqrt(u/xx);
     695         700 :       x = w * (a + c / (48*a*u) * (-5/u-c * (-10/xx + 6)));
     696             :     }
     697        1008 :     z = dbltor(x);
     698             :   }
     699             :   else
     700             :   { /* generic, hope for the best */
     701           0 :     long a = 4 * n - (J? 1: 3);
     702             :     GEN b, m;
     703           0 :     b = gmul(mppi(prec), gmul2n(gaddgs(gmul2n(nu,1), a), -2));
     704           0 :     m = gmul2n(gsqr(nu),2);
     705           0 :     z = gsub(b, gdiv(gsubgs(m, 1), gmul2n(b, 3)));
     706             :   }
     707        1008 :   return gc_GEN(av, besselrefine(z, nu, B, bit));
     708             : }
     709             : GEN
     710         511 : besseljzero(GEN nu, long k, long b) { return besselzero(nu, k, jbessel, b); }
     711             : GEN
     712         511 : besselyzero(GEN nu, long k, long b) { return besselzero(nu, k, ybessel, b); }
     713             : 
     714             : /***********************************************************************/
     715             : /**                    INCOMPLETE GAMMA FUNCTION                      **/
     716             : /***********************************************************************/
     717             : /* mx ~ |x|, b = bit accuracy */
     718             : static int
     719       16765 : gamma_use_asymp(GEN x, long b)
     720             : {
     721             :   long e;
     722       16765 :   if (is_real_t(typ(x)))
     723             :   {
     724       12733 :     pari_sp av = avma;
     725       12733 :     return gc_int(av, gcmpgs(R_abs_shallow(x), 3*b / 4) >= 0);
     726             :   }
     727        4032 :   e = gexpo(x); return e >= b || dblmodulus(x) >= 3*b / 4;
     728             : }
     729             : /* x a t_REAL */
     730             : static GEN
     731          28 : eint1r_asymp(GEN x, GEN expx, long prec)
     732             : {
     733          28 :   pari_sp av = avma, av2;
     734             :   GEN S, q, z, ix;
     735          28 :   long oldeq = LONG_MAX, esx = -prec2nbits(prec), j;
     736             : 
     737          28 :   if (realprec(x) < prec + EXTRAPREC64) x = rtor(x, prec+EXTRAPREC64);
     738          28 :   ix = invr(x); q = z = negr(ix);
     739          28 :   av2 = avma; S = addrs(q, 1);
     740          28 :   for (j = 2;; j++)
     741        1211 :   {
     742        1239 :     long eq = expo(q); if (eq < esx) break;
     743        1211 :     if ((j & 3) == 0)
     744             :     { /* guard against divergence */
     745         294 :       if (eq > oldeq) return gc_NULL(av); /* regressing, abort */
     746         294 :       oldeq = eq;
     747             :     }
     748        1211 :     q = mulrr(q, mulru(z, j)); S = addrr(S, q);
     749        1211 :     if (gc_needed(av2, 1)) (void)gc_all(av2, 2, &S, &q);
     750             :   }
     751          28 :   if (DEBUGLEVEL > 2) err_printf("eint1: using asymp\n");
     752          28 :   S = expx? divrr(S, expx): mulrr(S, mpexp(negr(x)));
     753          28 :   return gc_leaf(av, mulrr(S, ix));
     754             : }
     755             : /* cf incgam_asymp(0, x); z = -1/x
     756             :  *   exp(-x)/x * (1 + z + 2! z^2 + ...) */
     757             : static GEN
     758         105 : eint1_asymp(GEN x, GEN expx, long prec)
     759             : {
     760         105 :   pari_sp av = avma, av2;
     761             :   GEN S, q, z, ix;
     762         105 :   long oldeq = LONG_MAX, esx = -prec2nbits(prec), j;
     763             : 
     764         105 :   if (typ(x) != t_REAL) x = gtofp(x, prec+EXTRAPREC64);
     765         105 :   if (typ(x) == t_REAL) return eint1r_asymp(x, expx, prec);
     766         105 :   ix = ginv(x); q = z = gneg_i(ix);
     767         105 :   av2 = avma; S = gaddgs(q, 1);
     768         105 :   for (j = 2;; j++)
     769        5824 :   {
     770        5929 :     long eq = gexpo(q); if (eq < esx) break;
     771        5824 :     if ((j & 3) == 0)
     772             :     { /* guard against divergence */
     773        1442 :       if (eq > oldeq) return gc_NULL(av); /* regressing, abort */
     774        1442 :       oldeq = eq;
     775             :     }
     776        5824 :     q = gmul(q, gmulgu(z, j)); S = gadd(S, q);
     777        5824 :     if (gc_needed(av2, 1)) (void)gc_all(av2, 2, &S, &q);
     778             :   }
     779         105 :   if (DEBUGLEVEL > 2) err_printf("eint1: using asymp\n");
     780         105 :   S = expx? gdiv(S, expx): gmul(S, gexp(gneg_i(x), prec));
     781         105 :   return gc_upto(av, gmul(S, ix));
     782             : }
     783             : 
     784             : /* eint1(x) = incgam(0, x); typ(x) = t_REAL, x > 0 */
     785             : static GEN
     786        6524 : eint1p(GEN x, GEN expx)
     787             : {
     788             :   pari_sp av;
     789        6524 :   long prec = realprec(x), bit = prec2nbits(prec), i;
     790             :   double mx;
     791             :   GEN z, S, t, H, run;
     792             : 
     793        6524 :   if (gamma_use_asymp(x, bit)
     794          28 :       && (z = eint1r_asymp(x, expx, prec))) return z;
     795        6496 :   mx = rtodbl(x);
     796        6496 :   if (mx > 1)
     797        3591 :     prec += nbits2extraprec((mx+log(mx))/M_LN2 + 10);
     798             :   else
     799        2905 :     prec += EXTRAPREC64;
     800        6496 :   bit = prec2nbits(prec);
     801        6496 :   run = real_1(prec); x = rtor(x, prec);
     802        6496 :   av = avma; S = z = t = H = run;
     803      618178 :   for (i = 2; expo(S) - expo(t) <= bit; i++)
     804             :   {
     805      611682 :     H = addrr(H, divru(run,i)); /* H = sum_{k<=i} 1/k */
     806      611682 :     z = divru(mulrr(x,z), i);   /* z = x^(i-1)/i! */
     807      611682 :     t = mulrr(z, H); S = addrr(S, t);
     808      611682 :     if ((i & 0x1ff) == 0) (void)gc_all(av, 4, &z,&t,&S,&H);
     809             :   }
     810        6496 :   return subrr(mulrr(x, divrr(S,expx? expx: mpexp(x))),
     811             :                addrr(mplog(x), mpeuler(prec)));
     812             : }
     813             : /* eint1(x) = incgam(0, x); typ(x) = t_REAL, x < 0
     814             :  * rewritten from code contributed by Manfred Radimersky */
     815             : static GEN
     816         140 : eint1m(GEN x, GEN expx)
     817             : {
     818         140 :   GEN p1, q, S, y, z = cgetg(3, t_COMPLEX);
     819         140 :   long l  = realprec(x), n  = prec2nbits(l), j;
     820         140 :   pari_sp av = avma;
     821             : 
     822         140 :   y  = rtor(x, l + EXTRAPREC64); setsigne(y,1); /* |x| */
     823         140 :   if (gamma_use_asymp(y, n))
     824             :   { /* ~eint1_asymp: asymptotic expansion */
     825          14 :     p1 = q = invr(y); S = addrs(q, 1);
     826         560 :     for (j = 2; expo(q) >= -n; j++) {
     827         546 :       q = mulrr(q, mulru(p1, j));
     828         546 :       S = addrr(S, q);
     829             :     }
     830          14 :     y  = mulrr(p1, expx? divrr(S, expx): mulrr(S, mpexp(y)));
     831             :   }
     832             :   else
     833             :   {
     834         126 :     p1 = q = S = y;
     835       24248 :     for (j = 2; expo(q) - expo(S) >= -n; j++) {
     836       24122 :       p1 = mulrr(y, divru(p1, j)); /* (-x)^j/j! */
     837       24122 :       q = divru(p1, j);
     838       24122 :       S = addrr(S, q);
     839             :     }
     840         126 :     y  = addrr(S, addrr(logr_abs(x), mpeuler(l)));
     841             :   }
     842         140 :   y = gc_leaf(av, y); togglesign(y);
     843         140 :   gel(z, 1) = y;
     844         140 :   y = mppi(l); setsigne(y, -1);
     845         140 :   gel(z, 2) = y; return z;
     846             : }
     847             : 
     848             : /* real(z*log(z)-z), z = x+iy */
     849             : static double
     850        8372 : mygamma(double x, double y)
     851             : {
     852        8372 :   if (x == 0.) return -(M_PI/2)*fabs(y);
     853        8372 :   return (x/2)*log(x*x+y*y)-x-y*atan(y/x);
     854             : }
     855             : 
     856             : /* x^s exp(-x) */
     857             : static GEN
     858       10843 : expmx_xs(GEN s, GEN x, GEN logx, long prec)
     859             : {
     860             :   GEN z;
     861       10843 :   long ts = typ(s);
     862       10843 :   if (ts == t_INT || (ts == t_FRAC && absequaliu(gel(s,2), 2)))
     863        5264 :     z = gmul(gexp(gneg(x), prec), gpow(x, s, prec));
     864             :   else
     865        5579 :     z = gexp(gsub(gmul(s, logx? logx: glog(x,prec+EXTRAPREC64)), x), prec);
     866       10843 :   return z;
     867             : }
     868             : 
     869             : /* Not yet: doesn't work at low accuracy
     870             : #define INCGAM_CF
     871             : */
     872             : 
     873             : #ifdef INCGAM_CF
     874             : /* Is s very close to a nonpositive integer ? */
     875             : static int
     876             : isgammapole(GEN s, long bitprec)
     877             : {
     878             :   pari_sp av = avma;
     879             :   GEN t = imag_i(s);
     880             :   long e, b = bitprec - 10;
     881             : 
     882             :   if (gexpo(t) > - b) return 0;
     883             :   s = real_i(s);
     884             :   if (gsigne(s) > 0 && gexpo(s) > -b) return 0;
     885             :   (void)grndtoi(s, &e); return gc_bool(av, e < -b);
     886             : }
     887             : 
     888             : /* incgam using the continued fraction. x a t_REAL or t_COMPLEX, mx ~ |x|.
     889             :  * Assume precision(s), precision(x) >= prec */
     890             : static GEN
     891             : incgam_cf(GEN s, GEN x, double mx, long prec)
     892             : {
     893             :   GEN ms, y, S;
     894             :   long n, i, j, LS, bitprec = prec2nbits(prec);
     895             :   double rs, is, m;
     896             : 
     897             :   if (typ(s) == t_COMPLEX)
     898             :   {
     899             :     rs = gtodouble(gel(s,1));
     900             :     is = gtodouble(gel(s,2));
     901             :   }
     902             :   else
     903             :   {
     904             :     rs = gtodouble(s);
     905             :     is = 0.;
     906             :   }
     907             :   if (isgammapole(s, bitprec)) LS = 0;
     908             :   else
     909             :   {
     910             :     double bit,  LGS = mygamma(rs,is);
     911             :     LS = LGS <= 0 ? 0: ceil(LGS);
     912             :     bit = (LGS - (rs-1)*log(mx) + mx)/M_LN2;
     913             :     if (bit > 0)
     914             :     {
     915             :       prec += nbits2extraprec((long)bit);
     916             :       x = gtofp(x, prec);
     917             :       if (isinexactreal(s)) s = gtofp(s, prec);
     918             :     }
     919             :   }
     920             :   /* |ln(2*gamma(s)*sin(s*Pi))| <= ln(2) + |lngamma(s)| + |Im(s)*Pi|*/
     921             :   m = bitprec*M_LN2 + LS + M_LN2 + fabs(is)*M_PI + mx;
     922             :   if (rs < 1) m += (1 - rs)*log(mx);
     923             :   m /= 4;
     924             :   n = (long)(1 + m*m/mx);
     925             :   y = expmx_xs(gsubgs(s,1), x, NULL, prec);
     926             :   if (rs >= 0 && bitprec >= 512)
     927             :   {
     928             :     GEN A = cgetg(n+1, t_VEC), B = cgetg(n+1, t_VEC);
     929             :     ms = gsubsg(1, s);
     930             :     for (j = 1; j <= n; ++j)
     931             :     {
     932             :       gel(A,j) = ms;
     933             :       gel(B,j) = gmulsg(j, gsubgs(s,j));
     934             :       ms = gaddgs(ms, 2);
     935             :     }
     936             :     S = contfraceval_inv(mkvec2(A,B), x, -1);
     937             :   }
     938             :   else
     939             :   {
     940             :     GEN x_s = gsub(x, s);
     941             :     pari_sp av2 = avma;
     942             :     S = gdiv(gsubgs(s,n), gaddgs(x_s,n<<1));
     943             :     for (i=n-1; i >= 1; i--)
     944             :     {
     945             :       S = gdiv(gsubgs(s,i), gadd(gaddgs(x_s,i<<1),gmulsg(i,S)));
     946             :       if (gc_needed(av2,3))
     947             :       {
     948             :         if(DEBUGMEM>1) pari_warn(warnmem,"incgam_cf");
     949             :         S = gc_upto(av2, S);
     950             :       }
     951             :     }
     952             :     S = gaddgs(S,1);
     953             :   }
     954             :   return gmul(y, S);
     955             : }
     956             : #endif
     957             : 
     958             : static double
     959        6419 : findextraincgam(GEN s, GEN x)
     960             : {
     961        6419 :   double sig = gtodouble(real_i(s)), t = gtodouble(imag_i(s));
     962        6419 :   double xr = gtodouble(real_i(x)), xi = gtodouble(imag_i(x));
     963        6419 :   double exd = 0., Nx = xr*xr + xi*xi, D = Nx - t*t;
     964             :   long n;
     965             : 
     966        6419 :   if (xr < 0)
     967             :   {
     968         833 :     long ex = gexpo(x);
     969         833 :     if (ex > 0 && ex > gexpo(s)) exd = sqrt(Nx)*log(Nx)/2; /* |x| log |x| */
     970             :   }
     971        6419 :   if (D <= 0.) return exd;
     972        4977 :   n = (long)(sqrt(D)-sig);
     973        4977 :   if (n <= 0) return exd;
     974        1841 :   return maxdd(exd, (n*log(Nx)/2 - mygamma(sig+n, t) + mygamma(sig, t)) / M_LN2);
     975             : }
     976             : 
     977             : /* use exp(-x) * (x^s/s) * sum_{k >= 0} x^k / prod(i=1, k, s+i) */
     978             : static GEN
     979        6426 : incgamc_i(GEN s, GEN x, long *ptexd, long prec)
     980             : {
     981             :   GEN S, t, y;
     982             :   long l, n, i, exd;
     983        6426 :   pari_sp av = avma, av2;
     984             : 
     985        6426 :   if (gequal0(x))
     986             :   {
     987           7 :     if (ptexd) *ptexd = 0.;
     988           7 :     return gtofp(x, prec);
     989             :   }
     990        6419 :   l = precision(x);
     991        6419 :   if (!l) l = prec;
     992        6419 :   n = -prec2nbits(l)-1;
     993        6419 :   exd = (long)findextraincgam(s, x);
     994        6419 :   if (ptexd) *ptexd = exd;
     995        6419 :   if (exd > 0)
     996             :   {
     997        1666 :     long p = l + nbits2extraprec(exd);
     998        1666 :     x = gtofp(x, p);
     999        1666 :     if (isinexactreal(s)) s = gtofp(s, p);
    1000             :   }
    1001        4753 :   else x = gtofp(x, l+EXTRAPREC64);
    1002        6419 :   av2 = avma;
    1003        6419 :   S = gdiv(x, gaddsg(1,s));
    1004        6419 :   t = gaddsg(1, S);
    1005      770875 :   for (i=2; gexpo(S) >= n; i++)
    1006             :   {
    1007      764456 :     S = gdiv(gmul(x,S), gaddsg(i,s)); /* x^i / ((s+1)...(s+i)) */
    1008      764456 :     t = gadd(S,t);
    1009      764456 :     if (gc_needed(av2,3))
    1010             :     {
    1011           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"incgamc");
    1012           0 :       (void)gc_all(av2, 2, &S, &t);
    1013             :     }
    1014             :   }
    1015        6419 :   y = expmx_xs(s, x, NULL, prec);
    1016        6419 :   return gc_upto(av, gmul(gdiv(y,s), t));
    1017             : }
    1018             : 
    1019             : GEN
    1020        2226 : incgamc(GEN s, GEN x, long prec)
    1021        2226 : { return incgamc_i(s, x, NULL, prec); }
    1022             : 
    1023             : /* incgamma using asymptotic expansion:
    1024             :  *   exp(-x)x^(s-1)(1 + (s-1)/x + (s-1)(s-2)/x^2 + ...) */
    1025             : static GEN
    1026        2716 : incgam_asymp(GEN s, GEN x, long prec)
    1027             : {
    1028        2716 :   pari_sp av = avma, av2;
    1029             :   GEN S, q, cox, invx;
    1030        2716 :   long oldeq = LONG_MAX, eq, esx, j;
    1031        2716 :   int flint = (typ(s) == t_INT && signe(s) > 0);
    1032             : 
    1033        2716 :   x = gtofp(x,prec+EXTRAPREC64);
    1034        2716 :   invx = ginv(x);
    1035        2716 :   esx = -prec2nbits(prec);
    1036        2716 :   av2 = avma;
    1037        2716 :   q = gmul(gsubgs(s,1), invx);
    1038        2716 :   S = gaddgs(q, 1);
    1039        2716 :   for (j = 2;; j++)
    1040             :   {
    1041      123879 :     eq = gexpo(q); if (eq < esx) break;
    1042      121282 :     if (!flint && (j & 3) == 0)
    1043             :     { /* guard against divergence */
    1044       15778 :       if (eq > oldeq) return gc_NULL(av); /* regressing, abort */
    1045       15659 :       oldeq = eq;
    1046             :     }
    1047      121163 :     q = gmul(q, gmul(gsubgs(s,j), invx));
    1048      121163 :     S = gadd(S, q);
    1049      121163 :     if (gc_needed(av2, 1)) (void)gc_all(av2, 2, &S, &q);
    1050             :   }
    1051        2597 :   if (DEBUGLEVEL > 2) err_printf("incgam: using asymp\n");
    1052        2597 :   cox = expmx_xs(gsubgs(s,1), x, NULL, prec);
    1053        2597 :   return gc_upto(av, gmul(cox, S));
    1054             : }
    1055             : 
    1056             : /* gasx = incgam(s-n,x). Compute incgam(s,x)
    1057             :  * = (s-1)(s-2)...(s-n)gasx + exp(-x)x^(s-1) *
    1058             :  *   (1 + (s-1)/x + ... + (s-1)(s-2)...(s-n+1)/x^(n-1)) */
    1059             : static GEN
    1060         546 : incgam_asymp_partial(GEN s, GEN x, GEN gasx, long n, long prec)
    1061             : {
    1062             :   pari_sp av;
    1063         546 :   GEN S, q, cox, invx, s1 = gsubgs(s, 1), sprod;
    1064             :   long j;
    1065         546 :   cox = expmx_xs(s1, x, NULL, prec);
    1066         546 :   if (n == 1) return gadd(cox, gmul(s1, gasx));
    1067         546 :   invx = ginv(x);
    1068         546 :   av = avma;
    1069         546 :   q = gmul(s1, invx);
    1070         546 :   S = gaddgs(q, 1);
    1071       52164 :   for (j = 2; j < n; j++)
    1072             :   {
    1073       51618 :     q = gmul(q, gmul(gsubgs(s, j), invx));
    1074       51618 :     S = gadd(S, q);
    1075       51618 :     if (gc_needed(av, 2)) (void)gc_all(av, 2, &S, &q);
    1076             :   }
    1077         546 :   sprod = gmul(gmul(q, gpowgs(x, n-1)), gsubgs(s, n));
    1078         546 :   return gadd(gmul(cox, S), gmul(sprod, gasx));
    1079             : }
    1080             : 
    1081             : /* Assume s != 0; called when Re(s) <= 1/2 */
    1082             : static GEN
    1083        2401 : incgamspec(GEN s, GEN x, GEN g, long prec)
    1084             : {
    1085        2401 :   GEN q, S, cox = gen_0, P, sk, S1, S2, S3, F3, logx, mx;
    1086        2401 :   long n, esk, E, k = itos(ground(gneg(real_i(s)))); /* >= 0 */
    1087             : 
    1088        2401 :   if (k && gexpo(x) > 0)
    1089             :   {
    1090         245 :     GEN xk = gdivgu(x, k);
    1091         245 :     long bitprec = prec2nbits(prec);
    1092         245 :     double d = (gexpo(xk) > bitprec)? bitprec*M_LN2: log(dblmodulus(xk));
    1093         245 :     d = k * (d + 1) / M_LN2;
    1094         245 :     if (d > 0) prec += nbits2extraprec((long)d);
    1095         245 :     if (isinexactreal(s)) s = gtofp(s, prec);
    1096             :   }
    1097        2401 :   x = gtofp(x, maxss(precision(x), prec) + EXTRAPREC64);
    1098        2401 :   sk = gaddgs(s, k); /* |Re(sk)| <= 1/2 */
    1099        2401 :   logx = glog(x, prec);
    1100        2401 :   mx = gneg(x);
    1101        2401 :   if (k == 0) { S = gen_0; P = gen_1; }
    1102             :   else
    1103             :   {
    1104             :     long j;
    1105         854 :     q = ginv(s); S = q; P = s;
    1106       16926 :     for (j = 1; j < k; j++)
    1107             :     {
    1108       16072 :       GEN sj = gaddgs(s, j);
    1109       16072 :       q = gmul(q, gdiv(x, sj));
    1110       16072 :       S = gadd(S, q);
    1111       16072 :       P = gmul(P, sj);
    1112             :     }
    1113         854 :     cox = expmx_xs(s, x, logx, prec); /* x^s exp(-x) */
    1114         854 :     S = gmul(S, gneg(cox));
    1115             :   }
    1116        2401 :   if (k && gequal0(sk))
    1117         175 :     return gadd(S, gdiv(eint1(x, prec), P));
    1118        2226 :   esk = gexpo(sk);
    1119        2226 :   if (esk > -7)
    1120             :   {
    1121        1015 :     GEN a, b, PG = gmul(sk, P);
    1122        1015 :     if (g) g = gmul(g, PG);
    1123        1015 :     a = incgam0(gaddgs(sk,1), x, g, prec);
    1124        1015 :     if (k == 0) cox = expmx_xs(s, x, logx, prec);
    1125        1015 :     b = gmul(gpowgs(x, k), cox);
    1126        1015 :     return gadd(S, gdiv(gsub(a, b), PG));
    1127             :   }
    1128        1211 :   E = prec2nbits(prec) + 1;
    1129        1211 :   if (gexpo(x) > 0)
    1130             :   {
    1131         420 :     long X = (long)(dblmodulus(x)/M_LN2);
    1132         420 :     prec += 2*nbits2extraprec(X);
    1133         420 :     x = gtofp(x, prec); mx = gneg(x);
    1134         420 :     logx = glog(x, prec); sk = gtofp(sk, prec);
    1135         420 :     E += X;
    1136             :   }
    1137        1211 :   if (isinexactreal(sk)) sk = gtofp(sk, prec+EXTRAPREC64);
    1138             :   /* |sk| < 2^-7 is small, guard against cancellation */
    1139        1211 :   F3 = gexpm1(gmul(sk, logx), prec);
    1140             :   /* ( gamma(1+sk) - exp(sk log(x))) ) / sk */
    1141        1211 :   S1 = gdiv(gsub(ggamma1m1(sk, prec+EXTRAPREC64), F3), sk);
    1142        1211 :   q = x; S3 = gdiv(x, gaddsg(1,sk));
    1143      255523 :   for (n = 2; gexpo(q) - gexpo(S3) > -E; n++)
    1144             :   {
    1145      254312 :     q = gmul(q, gdivgu(mx, n));
    1146      254312 :     S3 = gadd(S3, gdiv(q, gaddsg(n, sk)));
    1147             :   }
    1148        1211 :   S2 = gadd(gadd(S1, S3), gmul(F3, S3));
    1149        1211 :   return gadd(S, gdiv(S2, P));
    1150             : }
    1151             : 
    1152             : /* return |x| */
    1153             : double
    1154    11826595 : dblmodulus(GEN x)
    1155             : {
    1156    11826595 :   if (typ(x) == t_COMPLEX)
    1157             :   {
    1158     1803571 :     double a = gtodouble(gel(x,1));
    1159     1803571 :     double b = gtodouble(gel(x,2));
    1160     1803571 :     return sqrt(a*a + b*b);
    1161             :   }
    1162             :   else
    1163    10023024 :     return fabs(gtodouble(x));
    1164             : }
    1165             : 
    1166             : /* Driver routine. If g != NULL, assume that g=gamma(s,prec). */
    1167             : GEN
    1168       11564 : incgam0(GEN s, GEN x, GEN g, long prec)
    1169             : {
    1170             :   pari_sp av;
    1171             :   long E, l;
    1172             :   GEN z, rs, is;
    1173             : 
    1174       11564 :   if (gequal0(x)) return g? gcopy(g): ggamma(s,prec);
    1175       11564 :   if (gequal0(s)) return eint1(x, prec);
    1176        9744 :   l = precision(s); if (!l) l = prec;
    1177        9744 :   E = prec2nbits(l);
    1178        9744 :   if (gamma_use_asymp(x, E) ||
    1179        8323 :       (typ(s) == t_INT && signe(s) > 0 && gexpo(x) >= expi(s)))
    1180        2716 :     if ((z = incgam_asymp(s, x, l))) return z;
    1181        7147 :   av = avma; E++;
    1182        7147 :   rs = real_i(s);
    1183        7147 :   is = imag_i(s);
    1184             : #ifdef INCGAM_CF
    1185             :   /* Can one use continued fraction ? */
    1186             :   if (gequal0(is) && gequal0(imag_i(x)) && gsigne(x) > 0)
    1187             :   {
    1188             :     double sd = gtodouble(rs), LB, UB;
    1189             :     double xd = gtodouble(real_i(x));
    1190             :     if (sd > 0) {
    1191             :       LB = 15 + 0.1205*E;
    1192             :       UB = 5 + 0.752*E;
    1193             :     } else {
    1194             :       LB = -6 + 0.1205*E;
    1195             :       UB = 5 + 0.752*E + fabs(sd)/54.;
    1196             :     }
    1197             :     if (xd >= LB && xd <= UB)
    1198             :     {
    1199             :       if (DEBUGLEVEL > 2) err_printf("incgam: using continued fraction\n");
    1200             :       return gc_upto(av, incgam_cf(s, x, xd, prec));
    1201             :     }
    1202             :   }
    1203             : #endif
    1204        7147 :   if (gsigne(rs) > 0 && gexpo(rs) >= -1)
    1205             :   { /* use complementary incomplete gamma */
    1206        4746 :     long n, egs, exd, precg, es = gexpo(s);
    1207        4746 :     if (es < 0) {
    1208         602 :       l += nbits2extraprec(-es) + 1;
    1209         602 :       x = gtofp(x, l);
    1210         602 :       if (isinexactreal(s)) s = gtofp(s, l);
    1211             :     }
    1212        4746 :     n = itos(gceil(rs));
    1213        4746 :     if (n > 100)
    1214             :     {
    1215             :       GEN gasx;
    1216         546 :       n -= 100;
    1217         546 :       if (es > 0)
    1218             :       {
    1219         546 :         es = mygamma(gtodouble(rs) - n, gtodouble(is)) / M_LN2;
    1220         546 :         if (es > 0)
    1221             :         {
    1222         546 :           l += nbits2extraprec(es);
    1223         546 :           x = gtofp(x, l);
    1224         546 :           if (isinexactreal(s)) s = gtofp(s, l);
    1225             :         }
    1226             :       }
    1227         546 :       gasx = incgam0(gsubgs(s, n), x, NULL, prec);
    1228         546 :       return gc_upto(av, incgam_asymp_partial(s, x, gasx, n, prec));
    1229             :     }
    1230        4200 :     if (DEBUGLEVEL > 2) err_printf("incgam: using power series 1\n");
    1231             :     /* egs ~ expo(gamma(s)) */
    1232        4200 :     precg = g? precision(g): 0;
    1233        4200 :     egs = g? gexpo(g): (long)(mygamma(gtodouble(rs), gtodouble(is)) / M_LN2);
    1234        4200 :     if (egs > 0) {
    1235        1946 :       l += nbits2extraprec(egs) + 1;
    1236        1946 :       x = gtofp(x, l);
    1237        1946 :       if (isinexactreal(s)) s = gtofp(s, l);
    1238        1946 :       if (precg < l) g = NULL;
    1239             :     }
    1240        4200 :     z = incgamc_i(s, x, &exd, l);
    1241        4200 :     if (exd > 0)
    1242             :     {
    1243         896 :       l += nbits2extraprec(exd);
    1244         896 :       if (isinexactreal(s)) s = gtofp(s, l);
    1245         896 :       if (precg < l) g = NULL;
    1246             :     }
    1247             :     else
    1248             :     { /* gamma(s) negligible ? Compute to lower accuracy */
    1249        3304 :       long e = gexpo(z) - egs;
    1250        3304 :       if (e > 3)
    1251             :       {
    1252         420 :         E -= e;
    1253         420 :         if (E <= 0) g = gen_0; else if (!g) g = ggamma(s, nbits2prec(E));
    1254             :       }
    1255             :     }
    1256             :     /* worry about possible cancellation */
    1257        4200 :     if (!g) g = ggamma(s, maxss(l,precision(z)));
    1258        4200 :     return gc_upto(av, gsub(g,z));
    1259             :   }
    1260        2401 :   if (DEBUGLEVEL > 2) err_printf("incgam: using power series 2\n");
    1261        2401 :   return gc_upto(av, incgamspec(s, x, g, l));
    1262             : }
    1263             : 
    1264             : GEN
    1265        1106 : incgam(GEN s, GEN x, long prec) { return incgam0(s, x, NULL, prec); }
    1266             : 
    1267             : /* x a t_REAL */
    1268             : GEN
    1269        2940 : mpeint1(GEN x, GEN expx)
    1270             : {
    1271        2940 :   long s = signe(x);
    1272             :   pari_sp av;
    1273             :   GEN z;
    1274        2940 :   if (!s) pari_err_DOMAIN("eint1", "x","=",gen_0, x);
    1275        2933 :   if (s < 0) return eint1m(x, expx);
    1276        2793 :   z = cgetr(realprec(x));
    1277        2793 :   av = avma; affrr(eint1p(x, expx), z);
    1278        2793 :   set_avma(av); return z;
    1279             : }
    1280             : 
    1281             : static GEN
    1282         357 : cxeint1(GEN x, long prec)
    1283             : {
    1284         357 :   pari_sp av = avma, av2;
    1285             :   GEN q, S, run, z, H;
    1286         357 :   long n, E = prec2nbits(prec);
    1287             : 
    1288         357 :   if (gamma_use_asymp(x, E) && (z = eint1_asymp(x, NULL, prec))) return z;
    1289         252 :   E++;
    1290         252 :   if (gexpo(x) > 0)
    1291             :   { /* take cancellation into account, log2(\sum |x|^n / n!) = |x| / log(2) */
    1292          42 :     double dbx = dblmodulus(x);
    1293          42 :     long X = (long)((dbx + log(dbx))/M_LN2 + 10);
    1294          42 :     prec += nbits2extraprec(X);
    1295          42 :     x = gtofp(x, prec); E += X;
    1296             :   }
    1297         252 :   if (DEBUGLEVEL > 2) err_printf("eint1: using power series\n");
    1298         252 :   run = real_1(prec);
    1299         252 :   av2 = avma;
    1300         252 :   S = z = q = H = run;
    1301       48384 :   for (n = 2; gexpo(q) - gexpo(S) >= -E; n++)
    1302             :   {
    1303       48132 :     H = addrr(H, divru(run, n)); /* H = sum_{k<=n} 1/k */
    1304       48132 :     z = gdivgu(gmul(x,z), n);   /* z = x^(n-1)/n! */
    1305       48132 :     q = gmul(z, H); S = gadd(S, q);
    1306       48132 :     if ((n & 0x1ff) == 0) (void)gc_all(av2, 4, &z, &q, &S, &H);
    1307             :   }
    1308         252 :   S = gmul(gmul(x, S), gexp(gneg_i(x), prec));
    1309         252 :   return gc_upto(av, gsub(S, gadd(glog(x, prec), mpeuler(prec))));
    1310             : }
    1311             : 
    1312             : GEN
    1313        3297 : eint1(GEN x, long prec)
    1314             : {
    1315        3297 :   switch(typ(x))
    1316             :   {
    1317         357 :     case t_COMPLEX: return cxeint1(x, prec);
    1318        2541 :     case t_REAL: break;
    1319         399 :     default: x = gtofp(x, prec);
    1320             :   }
    1321        2940 :   return mpeint1(x,NULL);
    1322             : }
    1323             : 
    1324             : GEN
    1325          49 : veceint1(GEN C, GEN nmax, long prec)
    1326             : {
    1327          49 :   if (!nmax) return eint1(C,prec);
    1328           7 :   if (typ(nmax) != t_INT) pari_err_TYPE("veceint1",nmax);
    1329           7 :   if (typ(C) != t_REAL) {
    1330           7 :     C = gtofp(C, prec);
    1331           7 :     if (typ(C) != t_REAL) pari_err_TYPE("veceint1",C);
    1332             :   }
    1333           7 :   if (signe(C) <= 0) pari_err_DOMAIN("veceint1", "argument", "<=", gen_0,C);
    1334           7 :   return mpveceint1(C, NULL, itos(nmax));
    1335             : }
    1336             : 
    1337             : /* j > 0, a t_REAL. Return sum_{m >= 0} a^m / j(j+1)...(j+m)).
    1338             :  * Stop when expo(summand) < E; note that s(j-1) = (a s(j) + 1) / (j-1). */
    1339             : static GEN
    1340         231 : mp_sum_j(GEN a, long j, long E, long prec)
    1341             : {
    1342         231 :   pari_sp av = avma;
    1343         231 :   GEN q = divru(real_1(prec), j), s = q;
    1344             :   long m;
    1345        4290 :   for (m = 0;; m++)
    1346             :   {
    1347        4290 :     if (expo(q) < E) break;
    1348        4059 :     q = mulrr(q, divru(a, m+j));
    1349        4059 :     s = addrr(s, q);
    1350             :   }
    1351         231 :   return gc_leaf(av, s);
    1352             : }
    1353             : /* Return the s_a(j), j <= J */
    1354             : static GEN
    1355         231 : sum_jall(GEN a, long J, long prec)
    1356             : {
    1357         231 :   GEN s = cgetg(J+1, t_VEC);
    1358         231 :   long j, E = -prec2nbits(prec) - 5;
    1359         231 :   gel(s, J) = mp_sum_j(a, J, E, prec);
    1360        9624 :   for (j = J-1; j; j--)
    1361        9393 :     gel(s,j) = divru(addrs(mulrr(a, gel(s,j+1)), 1), j);
    1362         231 :   return s;
    1363             : }
    1364             : 
    1365             : /* T a dense t_POL with t_REAL coeffs. Return T(n) [faster than poleval] */
    1366             : static GEN
    1367      364903 : rX_s_eval(GEN T, long n)
    1368             : {
    1369      364903 :   long i = lg(T)-1;
    1370      364903 :   GEN c = gel(T,i);
    1371     7536888 :   for (i--; i>=2; i--) c = gadd(mulrs(c,n),gel(T,i));
    1372      364903 :   return c;
    1373             : }
    1374             : 
    1375             : /* C>0 t_REAL, eC = exp(C). Return eint1(n*C) for 1<=n<=N. Absolute accuracy */
    1376             : GEN
    1377         238 : mpveceint1(GEN C, GEN eC, long N)
    1378             : {
    1379         238 :   const long prec = realprec(C);
    1380         238 :   long Nmin = 15; /* >= 1. E.g. between 10 and 30, but little effect */
    1381         238 :   GEN en, v, w = cgetg(N+1, t_VEC);
    1382             :   pari_sp av0;
    1383             :   double DL;
    1384             :   long n, j, jmax, jmin;
    1385         238 :   if (!N) return w;
    1386      368641 :   for (n = 1; n <= N; n++) gel(w,n) = cgetr(prec);
    1387         238 :   av0 = avma;
    1388         238 :   if (N < Nmin) Nmin = N;
    1389         238 :   if (!eC) eC = mpexp(C);
    1390         238 :   en = eC; affrr(eint1p(C, en), gel(w,1));
    1391        3500 :   for (n = 2; n <= Nmin; n++)
    1392             :   {
    1393             :     pari_sp av2;
    1394        3262 :     en = mulrr(en,eC); /* exp(n C) */
    1395        3262 :     av2 = avma;
    1396        3262 :     affrr(eint1p(mulru(C,n), en), gel(w,n));
    1397        3262 :     set_avma(av2);
    1398             :   }
    1399         238 :   if (Nmin == N) { set_avma(av0); return w; }
    1400             : 
    1401         231 :   DL = prec2nbits_mul(prec, M_LN2) + 5;
    1402         231 :   jmin = ceil(DL/log((double)N)) + 1;
    1403         231 :   jmax = ceil(DL/log((double)Nmin)) + 1;
    1404         231 :   v = sum_jall(C, jmax, prec);
    1405         231 :   en = powrs(eC, -N); /* exp(-N C) */
    1406         231 :   affrr(eint1p(mulru(C,N), invr(en)), gel(w,N));
    1407        6041 :   for (j = jmin, n = N-1; j <= jmax; j++)
    1408             :   {
    1409        5810 :     long limN = maxss((long)ceil(exp(DL/j)), Nmin);
    1410             :     GEN polsh;
    1411        5810 :     setlg(v, j+1);
    1412        5810 :     polsh = RgV_to_RgX_reverse(v, 0);
    1413      370713 :     for (; n >= limN; n--)
    1414             :     {
    1415      364903 :       pari_sp av2 = avma;
    1416      364903 :       GEN S = divri(mulrr(en, rX_s_eval(polsh, -n)), powuu(n,j));
    1417             :       /* w[n+1] - exp(-n C) * polsh(-n) / (-n)^j */
    1418      364903 :       GEN c = odd(j)? addrr(gel(w,n+1), S) : subrr(gel(w,n+1), S);
    1419      364903 :       affrr(c, gel(w,n)); set_avma(av2);
    1420      364903 :       en = mulrr(en,eC); /* exp(-n C) */
    1421             :     }
    1422             :   }
    1423         231 :   set_avma(av0); return w;
    1424             : }
    1425             : 
    1426             : /* erfc via numerical integration : assume real(x)>=1 */
    1427             : static GEN
    1428          14 : cxerfc_r1(GEN x, long prec)
    1429             : {
    1430             :   GEN h, h2, eh2, denom, res, lambda;
    1431             :   long u, v;
    1432          14 :   const double D = prec2nbits_mul(prec, M_LN2);
    1433          14 :   const long npoints = (long)ceil(D/M_PI)+1;
    1434          14 :   pari_sp av = avma;
    1435             :   {
    1436          14 :     double t = exp(-2*M_PI*M_PI/D); /* ~exp(-2*h^2) */
    1437          14 :     v = 30; /* bits that fit in both long and double mantissa */
    1438          14 :     u = (long)floor(t*(1L<<v));
    1439             :     /* define exp(-2*h^2) to be u*2^(-v) */
    1440             :   }
    1441          14 :   incrprec(prec);
    1442          14 :   x = gtofp(x,prec);
    1443          14 :   eh2 = sqrtr_abs(rtor(shiftr(dbltor(u),-v),prec));
    1444          14 :   h2 = negr(logr_abs(eh2));
    1445          14 :   h = sqrtr_abs(h2);
    1446          14 :   lambda = gdiv(x,h);
    1447          14 :   denom = gsqr(lambda);
    1448             :   { /* res = h/x + 2*x*h*sum(k=1,npoints,exp(-(k*h)^2)/(lambda^2+k^2)); */
    1449             :     GEN Uk; /* = exp(-(kh)^2) */
    1450          14 :     GEN Vk = eh2;/* = exp(-(2k+1)h^2) */
    1451          14 :     pari_sp av2 = avma;
    1452             :     long k;
    1453             :     /* k = 0 moved out for efficiency */
    1454          14 :     denom = gaddsg(1,denom);
    1455          14 :     Uk = Vk;
    1456          14 :     Vk = mulur(u,Vk); shiftr_inplace(Vk, -v);
    1457          14 :     res = gdiv(Uk, denom);
    1458         420 :     for (k = 1; k < npoints; k++)
    1459             :     {
    1460         406 :       if ((k & 255) == 0) (void)gc_all(av2,4,&denom,&Uk,&Vk,&res);
    1461         406 :       denom = gaddsg(2*k+1,denom);
    1462         406 :       Uk = mpmul(Uk,Vk);
    1463         406 :       Vk = mulur(u,Vk); shiftr_inplace(Vk, -v);
    1464         406 :       res = gadd(res, gdiv(Uk, denom));
    1465             :     }
    1466             :   }
    1467          14 :   res = gmul(res, gshift(lambda,1));
    1468             :   /* 0 term : */
    1469          14 :   res = gadd(res, ginv(lambda));
    1470          14 :   res = gmul(res, gdiv(gexp(gneg(gsqr(x)), prec), mppi(prec)));
    1471          14 :   if (rtodbl(real_i(x)) < sqrt(D))
    1472             :   {
    1473          14 :     GEN t = gmul(divrr(Pi2n(1,prec),h), x);
    1474          14 :     res = gsub(res, gdivsg(2, cxexpm1(t, prec)));
    1475             :   }
    1476          14 :   return gc_upto(av,res);
    1477             : }
    1478             : 
    1479             : static GEN
    1480           7 : sererfc(GEN x, long prec)
    1481             : {
    1482           7 :   GEN u, z = invr(sqrtr_abs(Pi2n(-2,prec)));
    1483           7 :   setsigne(z, -1); /* -2/sqrt(Pi) */
    1484           7 :   z = gmul(z, integser(gmul(derivser(x), gexp(gneg(gsqr(x)), prec))));
    1485           7 :   u = polcoef_i(x, 0, varn(x));
    1486           7 :   if (!gcmp0(u)) z = gadd(z, gerfc(u, prec));
    1487           7 :   return z;
    1488             : }
    1489             : 
    1490             : GEN
    1491          70 : gerfc(GEN x, long prec)
    1492             : {
    1493             :   GEN z, xr, xi, res;
    1494             :   long s;
    1495             :   pari_sp av;
    1496             : 
    1497          70 :   switch(typ(x))
    1498             :   {
    1499          63 :     case t_INT: case t_REAL: case t_FRAC: case t_COMPLEX:
    1500          63 :       break;
    1501           7 :     default:
    1502           7 :       av = avma;
    1503           7 :       if ((z = toser_i(x))) return gc_upto(av, sererfc(z,prec));
    1504           0 :       return trans_eval("erfc",gerfc,x,prec);
    1505             :   }
    1506             :   /* x a complex scalar */
    1507          63 :   x = trans_fix_arg(&prec,&x,&xr,&xi, &av,&res);
    1508          63 :   s = signe(xr);
    1509          63 :   if (s > 0 || (s == 0 && signe(xi) >= 0)) {
    1510          49 :     if (cmprs(xr, 1) > 0) /* use numerical integration */
    1511          14 :       z = cxerfc_r1(x, prec);
    1512             :     else
    1513             :     { /* erfc(x) = incgam(1/2,x^2)/sqrt(Pi) */
    1514          35 :       GEN sqrtpi = sqrtr(mppi(prec));
    1515          35 :       z = incgam0(ghalf, gsqr(x), sqrtpi, prec);
    1516          35 :       z = gdiv(z, sqrtpi);
    1517             :     }
    1518             :   }
    1519             :   else
    1520             :   { /* erfc(-x)=2-erfc(x) */
    1521             :     /* FIXME could decrease prec
    1522             :     long size = nbits2extraprec((imag(x)^2-real(x)^2)/log(2));
    1523             :     prec = size > 0 ? prec : prec + size;
    1524             :     */
    1525             :     /* NOT gsubsg(2, ...) : would create a result of
    1526             :      * huge accuracy if re(x)>>1, rounded to 2 by subsequent affc_fixlg... */
    1527          14 :     z = gsub(real2n(1,prec+EXTRAPREC64), gerfc(gneg(x), prec));
    1528             :   }
    1529          63 :   set_avma(av); return affc_fixlg(z, res);
    1530             : }
    1531             : 
    1532             : /***********************************************************************/
    1533             : /**                                                                   **/
    1534             : /**                      RIEMANN ZETA FUNCTION                        **/
    1535             : /**                                                                   **/
    1536             : /***********************************************************************/
    1537             : static const double log2PI = 1.83787706641;
    1538             : 
    1539             : static double
    1540        4585 : get_xinf(double beta)
    1541             : {
    1542        4585 :   const double maxbeta = 0.06415003; /* 3^(-2.5) */
    1543             :   double x0, y0, x1;
    1544             : 
    1545        4585 :   if (beta < maxbeta) return beta + pow(3*beta, 1.0/3.0);
    1546        4585 :   x0 = beta + M_PI/2.0;
    1547             :   for(;;)
    1548             :   {
    1549        7539 :     y0 = x0*x0;
    1550        7539 :     x1 = (beta+atan(x0)) * (1+y0) / y0 - 1/x0;
    1551        7539 :     if (0.99*x0 < x1) return x1;
    1552        2954 :     x0 = x1;
    1553             :   }
    1554             : }
    1555             : /* optimize for zeta( s + it, prec ), assume |s-1| > 0.1
    1556             :  * (if gexpo(u = s-1) < -5, we use the functional equation s->1-s) */
    1557             : static int
    1558       18494 : optim_zeta(GEN S, long prec, long *pp, long *pn)
    1559             : {
    1560             :   double s, t, alpha, beta, n, B;
    1561             :   long p;
    1562       18494 :   if (typ(S) == t_REAL) {
    1563       10129 :     t = 0.;
    1564       10129 :     s = rtodbl(S);
    1565             :   } else {
    1566        8365 :     t = fabs( rtodbl(gel(S,2)) );
    1567        8365 :     if (t > 2500) return 0; /* lfunlarge */
    1568        8316 :     s = rtodbl(gel(S,1));
    1569             :   }
    1570             : 
    1571       18445 :   B = prec2nbits_mul(prec, M_LN2);
    1572       18445 :   if (s > 0 && !t) /* positive real input */
    1573             :   {
    1574       10080 :     beta = B + 0.61 + s*(log2PI - log(s));
    1575       10080 :     if (beta > 0)
    1576             :     {
    1577        5110 :       p = (long)ceil(beta / 2.0);
    1578        5110 :       n = fabs(s + 2*p-1)/(2*M_PI);
    1579             :     }
    1580             :     else
    1581             :     {
    1582        4970 :       p = 0;
    1583        4970 :       n = exp((B - M_LN2) / s);
    1584             :     }
    1585             :   }
    1586        8365 :   else if (s <= 0 || t < 0.01) /* s < 0 may occur if s ~ 0 */
    1587        3773 :   { /* TODO: the crude bounds below are generally valid. Optimize ? */
    1588        3773 :     double l,l2, la = 1.; /* heuristic */
    1589        3773 :     double rlog, ilog; dblclog(s-1,t, &rlog,&ilog);
    1590        3773 :     l2 = (s - 0.5)*rlog - t*ilog; /* = Re( (S - 1/2) log (S-1) ) */
    1591        3773 :     l = (B - l2 + s*log2PI) / (2. * (1.+ log((double)la)));
    1592        3773 :     l2 = dblcabs(s, t)/2;
    1593        3773 :     if (l < l2) l = l2;
    1594        3773 :     p = (long) ceil(l); if (p < 2) p = 2;
    1595        3773 :     n = 1 + dblcabs(p+s/2.-.25, t/2) * la / M_PI;
    1596             :   }
    1597             :   else
    1598             :   {
    1599        4592 :     double sn = dblcabs(s, t), L = log(sn/s);
    1600        4592 :     alpha = B - 0.39 + L + s*(log2PI - log(sn));
    1601        4592 :     beta = (alpha+s)/t - atan(s/t);
    1602        4592 :     p = 0;
    1603        4592 :     if (beta > 0)
    1604             :     {
    1605        4585 :       beta = 1.0 - s + t * get_xinf(beta);
    1606        4585 :       if (beta > 0) p = (long)ceil(beta / 2.0);
    1607             :     }
    1608             :     else
    1609           7 :       if (s < 1.0) p = 1;
    1610        4592 :     n = p? dblcabs(s + 2*p-1, t) / (2*M_PI) : exp((B-M_LN2+L) / s);
    1611             :   }
    1612       18445 :   *pp = p;
    1613       18445 :   *pn = (long)ceil(n);
    1614       18445 :   if (*pp < 0 || *pn < 0) pari_err_OVERFLOW("zeta");
    1615       18445 :   return 1;
    1616             : }
    1617             : 
    1618             : /* zeta(a*j+b), j=0..N-1, b>1, using sumalt. Johansonn's thesis, Algo 4.7.1 */
    1619             : static GEN
    1620         705 : veczetas(long a, long b, long N, long prec)
    1621             : {
    1622         705 :   const long n = ceil(2 + prec2nbits_mul(prec, M_LN2/1.7627));
    1623         705 :   pari_sp av = avma;
    1624         705 :   GEN c, d, z = zerovec(N);
    1625             :   long j, k;
    1626         705 :   c = d = int2n(2*n-1);
    1627       88479 :   for (k = n; k > 1; k--)
    1628             :   {
    1629       87774 :     GEN u, t = divii(d, powuu(k,b));
    1630       87751 :     if (!odd(k)) t = negi(t);
    1631       87783 :     gel(z,1) = addii(gel(z,1), t);
    1632       87842 :     u = powuu(k,a);
    1633     4113900 :     for (j = 1; j < N; j++)
    1634             :     {
    1635     4070067 :       t = divii(t,u); if (!signe(t)) break;
    1636     4022233 :       gel(z,j+1) = addii(gel(z,j+1), t);
    1637             :     }
    1638       89202 :     c = muluui(k,2*k-1,c);
    1639       87792 :     c = diviuuexact(c, 2*(n-k+1),n+k-1);
    1640       87767 :     d = addii(d,c);
    1641       87775 :     if (gc_needed(av,3))
    1642             :     {
    1643           7 :       if(DEBUGMEM>1) pari_warn(warnmem,"zetaBorwein, k = %ld", k);
    1644           7 :       (void)gc_all(av, 3, &c,&d,&z);
    1645             :     }
    1646             :   }
    1647             :   /* k = 1 */
    1648       53410 :   for (j = 1; j <= N; j++) gel(z,j) = addii(gel(z,j), d);
    1649         707 :   d = addiu(d, 1);
    1650       53423 :   for (j = 0, k = b - 1; j < N; j++, k += a)
    1651       52718 :     gel(z,j+1) = rdivii(shifti(gel(z,j+1), k), subii(shifti(d,k), d), prec);
    1652         705 :   return z;
    1653             : }
    1654             : /* zeta(a*j+b), j=0..N-1, b > 1, a*(N-1) + b > 1, using sumalt.
    1655             :  * a <= 0 is allowed (including the silly a = 0) */
    1656             : GEN
    1657         224 : veczeta(GEN a, GEN b, long N, long prec)
    1658             : {
    1659         224 :   pari_sp av = avma;
    1660             :   long n, j, k;
    1661             :   GEN L, c, d, z;
    1662         224 :   if (typ(a) == t_INT && typ(b) == t_INT)
    1663         126 :     return gc_GEN(av, veczetas(itos(a),  itos(b), N, prec));
    1664          98 :   z = zerovec(N);
    1665          98 :   n = ceil(2 + prec2nbits_mul(prec, M_LN2/1.7627));
    1666          98 :   c = d = int2n(2*n-1);
    1667       14197 :   for (k = n; k; k--)
    1668             :   {
    1669             :     GEN u, t;
    1670       14099 :     L = logr_abs(utor(k, prec)); /* log(k) */
    1671       14099 :     t = gdiv(d, gexp(gmul(b, L), prec)); /* d / k^b */
    1672       14099 :     if (!odd(k)) t = gneg(t);
    1673       14099 :     gel(z,1) = gadd(gel(z,1), t);
    1674       14099 :     u = gexp(gmul(a, L), prec);
    1675      898001 :     for (j = 1; j < N; j++)
    1676             :     {
    1677      890325 :       t = gdiv(t,u); if (gexpo(t) < 0) break;
    1678      883902 :       gel(z,j+1) = gadd(gel(z,j+1), t);
    1679             :     }
    1680       14099 :     c = muluui(k,2*k-1,c);
    1681       14099 :     c = diviuuexact(c, 2*(n-k+1),n+k-1);
    1682       14099 :     d = addii(d,c);
    1683       14099 :     if (gc_needed(av,3))
    1684             :     {
    1685           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"veczeta, k = %ld", k);
    1686           0 :       (void)gc_all(av, 3, &c,&d,&z);
    1687             :     }
    1688             :   }
    1689          98 :   L = mplog2(prec);
    1690        5824 :   for (j = 0; j < N; j++)
    1691             :   {
    1692        5726 :     GEN u = gsubgs(gadd(b, gmulgu(a,j)), 1);
    1693        5726 :     GEN w = gexp(gmul(u, L), prec); /* 2^u */
    1694        5726 :     gel(z,j+1) = gdiv(gmul(gel(z,j+1), w), gmul(d,gsubgs(w,1)));
    1695             :   }
    1696          98 :   return gc_GEN(av, z);
    1697             : }
    1698             : 
    1699             : GEN
    1700       27443 : constzeta(long n, long prec)
    1701             : {
    1702       27443 :   GEN o = zetazone, z;
    1703       27443 :   long l = o? lg(o): 0;
    1704             :   pari_sp av;
    1705       27443 :   if (l > n)
    1706             :   {
    1707       27035 :     long p = realprec(gel(o,1));
    1708       27035 :     if (p >= prec) return o;
    1709             :   }
    1710         579 :   n = maxss(n, l + 15);
    1711         579 :   av = avma; z = veczetas(1, 2, n-1, prec);
    1712         579 :   zetazone = gclone(vec_prepend(z, mpeuler(prec)));
    1713         578 :   set_avma(av); guncloneNULL(o); return zetazone;
    1714             : }
    1715             : 
    1716             : /* zeta(s) using sumalt, case h=0,N=1. Assume s > 1 */
    1717             : static GEN
    1718         598 : zetaBorwein(long s, long prec)
    1719             : {
    1720         598 :   pari_sp av = avma;
    1721         598 :   const long n = ceil(2 + prec2nbits_mul(prec, M_LN2/1.7627));
    1722             :   long k;
    1723         598 :   GEN c, d, z = gen_0;
    1724         598 :   c = d = int2n(2*n-1);
    1725      135672 :   for (k = n; k; k--)
    1726             :   {
    1727      135074 :     GEN t = divii(d, powuu(k,s));
    1728      135074 :     z = odd(k)? addii(z,t): subii(z,t);
    1729      135074 :     c = muluui(k,2*k-1,c);
    1730      135074 :     c = diviuuexact(c, 2*(n-k+1),n+k-1);
    1731      135074 :     d = addii(d,c);
    1732      135074 :     if (gc_needed(av,3))
    1733             :     {
    1734           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"zetaBorwein, k = %ld", k);
    1735           0 :       (void)gc_all(av, 3, &c,&d,&z);
    1736             :     }
    1737             :   }
    1738         598 :   return rdivii(shifti(z, s-1), subii(shifti(d,s-1), d), prec);
    1739             : }
    1740             : 
    1741             : /* assume k != 1 */
    1742             : GEN
    1743        4760 : szeta(long k, long prec)
    1744             : {
    1745        4760 :   pari_sp av = avma;
    1746             :   GEN z;
    1747             : 
    1748        4760 :   if (!k) { z = real2n(-1, prec); setsigne(z,-1); return z; }
    1749        4753 :   if (k < 0)
    1750             :   {
    1751           0 :     if (!odd(k)) return gen_0;
    1752             :     /* the one value such that k < 0 and 1 - k < 0, due to overflow */
    1753           0 :     if ((ulong)k == (HIGHBIT | 1))
    1754           0 :       pari_err_OVERFLOW("zeta [large negative argument]");
    1755           0 :     k = 1-k;
    1756           0 :     z = bernreal(k, prec); togglesign(z);
    1757           0 :     return gc_leaf(av, divru(z, k));
    1758             :   }
    1759             :   /* k > 1 */
    1760        4753 :   if (k > prec2nbits(prec)+1) return real_1(prec);
    1761        4746 :   if (zetazone && realprec(gel(zetazone,1)) >= prec && lg(zetazone) > k)
    1762         140 :     return rtor(gel(zetazone, k), prec);
    1763        4606 :   if (!odd(k))
    1764             :   {
    1765             :     GEN B;
    1766        2835 :     if (!bernzone) constbern(0);
    1767        2835 :     if (k < lg(bernzone))
    1768        2380 :       B = gel(bernzone, k>>1);
    1769             :     else
    1770             :     {
    1771         455 :       if (bernbitprec(k) > prec2nbits(prec))
    1772           0 :         return gc_upto(av, invr(inv_szeta_euler(k, prec)));
    1773         455 :       B = bernfrac(k);
    1774             :     }
    1775             :     /* B = B_k */
    1776        2835 :     z = gmul(powru(Pi2n(1, prec + EXTRAPREC64), k), B);
    1777        2835 :     z = k < 410? rtor(divri(z, mpfact(k)), prec)
    1778        2835 :                : divrr(z, mpfactr(k,prec));
    1779        2835 :     setsigne(z, 1); shiftr_inplace(z, -1);
    1780             :   }
    1781             :   else
    1782             :   {
    1783        1771 :     double p = prec2nbits_mul(prec,0.393); /* bit / log_2(3+sqrt(8)) */
    1784        1771 :     p = log2(p * log(p));
    1785        3542 :     z = (p * k > prec2nbits(prec))? invr(inv_szeta_euler(k, prec))
    1786        1771 :                                   : zetaBorwein(k, prec);
    1787             :   }
    1788        4606 :   return gc_leaf(av, z);
    1789             : }
    1790             : 
    1791             : /* Ensure |1-s| >= 1/32 and (|s| <= 1/32 or real(s) >= 1/2) */
    1792             : static int
    1793       18508 : zeta_funeq(GEN *ps)
    1794             : {
    1795       18508 :   GEN s = *ps, u;
    1796       18508 :   if (typ(s) == t_REAL)
    1797             :   {
    1798       10136 :     u = subsr(1, s);
    1799       10136 :     if (expo(u) >= -5
    1800       10108 :         && ((signe(s) > 0 && expo(s) >= -1) || expo(s) <= -5)) return 0;
    1801             :   }
    1802             :   else
    1803             :   {
    1804        8372 :     GEN sig = gel(s,1);
    1805        8372 :     if (fabs(rtodbl(gel(s,2))) > 2500) return 0; /* lfunlarge */
    1806        8323 :     u = gsubsg(1, s);
    1807        8323 :     if (gexpo(u) >= -5
    1808        8316 :         && ((signe(sig) > 0 && expo(sig) >= -1) || gexpo(s) <= -5)) return 0;
    1809             :   }
    1810        1526 :   *ps = u; return 1;
    1811             : }
    1812             : /* s0 a t_INT, t_REAL or t_COMPLEX.
    1813             :  * If a t_INT, assume it's not a trivial case (i.e we have s0 > 1, odd) */
    1814             : static GEN
    1815       18515 : czeta(GEN s0, long prec)
    1816             : {
    1817       18515 :   GEN ms, s, u, y, res, tes, sig, tau, invn2, ns, Ns, funeq_factor = NULL;
    1818             :   long i, nn, lim, lim2;
    1819       18515 :   pari_sp av0 = avma, av, av2;
    1820             :   pari_timer T;
    1821             : 
    1822       18515 :   if (DEBUGLEVEL>2) timer_start(&T);
    1823       18515 :   s = trans_fix_arg(&prec,&s0,&sig,&tau,&av,&res);
    1824       18515 :   if (typ(s0) == t_INT) return gc_upto(av0, gzeta(s0, prec));
    1825       18508 :   if (zeta_funeq(&s)) /* s -> 1-s */
    1826             :   { /* Gamma(s) (2Pi)^-s 2 cos(Pi s/2) [new s] */
    1827        1526 :     GEN t = gmul(ggamma(s,prec), pow2Pis(gsubgs(s0,1), prec));
    1828        1526 :     sig = real_i(s);
    1829        1526 :     funeq_factor = gmul2n(gmul(t, gsin(gmul(Pi2n(-1,prec),s0), prec)), 1);
    1830             :   }
    1831       18508 :   if (gcmpgs(sig, prec2nbits(prec) + 1) > 0) { /* zeta(s) = 1 */
    1832          14 :     if (!funeq_factor) { set_avma(av0); return real_1(prec); }
    1833           7 :     return gc_upto(av0, funeq_factor);
    1834             :   }
    1835       18494 :   if (!optim_zeta(s, prec, &lim, &nn))
    1836             :   {
    1837          49 :     long bit = prec2nbits(prec);
    1838          49 :     y = lfun(lfuninit(gen_1, cgetg(1,t_VEC), 0, bit), s, bit);
    1839          49 :     if (funeq_factor) y = gmul(y, funeq_factor);
    1840          49 :     set_avma(av); return affc_fixlg(y,res);
    1841             :   }
    1842       18445 :   if (DEBUGLEVEL>2) err_printf("lim, nn: [%ld, %ld]\n", lim, nn);
    1843       18445 :   ms = gneg(s);
    1844       18445 :   if (umuluu_le(nn, prec, 10000000))
    1845             :   {
    1846       18445 :     incrprec(prec); /* one extra word of precision */
    1847       18445 :     Ns = vecpowug(nn, ms, prec);
    1848       18445 :     ns = gel(Ns,nn); setlg(Ns, nn);
    1849       18445 :     y = gadd(gmul2n(ns, -1), RgV_sum(Ns));
    1850             :   }
    1851             :   else
    1852             :   {
    1853           0 :     Ns = dirpowerssum(nn, ms, 0, prec);
    1854           0 :     incrprec(prec); /* one extra word of precision */
    1855           0 :     ns = gpow(utor(nn, prec), ms, prec);
    1856           0 :     y = gsub(Ns, gmul2n(ns, -1));
    1857             :   }
    1858       18445 :   if (DEBUGLEVEL>2) timer_printf(&T,"sum from 1 to N");
    1859       18445 :   constbern(lim);
    1860       18445 :   if (DEBUGLEVEL>2) timer_start(&T);
    1861       18445 :   invn2 = divri(real_1(prec), sqru(nn)); lim2 = lim<<1;
    1862       18445 :   tes = bernfrac(lim2);
    1863             :   {
    1864             :     GEN s1, s2, s3, s4, s5;
    1865       18445 :     s2 = gmul(s, gsubgs(s,1));
    1866       18445 :     s3 = gmul2n(invn2,3);
    1867       18445 :     av2 = avma;
    1868       18445 :     s1 = gsubgs(gmul2n(s,1), 1);
    1869       18445 :     s4 = gmul(invn2, gmul2n(gaddsg(4*lim-2,s1),1));
    1870       18445 :     s5 = gmul(invn2, gadd(s2, gmulsg(lim2, gaddgs(s1, lim2))));
    1871     1584174 :     for (i = lim2-2; i>=2; i -= 2)
    1872             :     {
    1873     1565729 :       s5 = gsub(s5, s4);
    1874     1565729 :       s4 = gsub(s4, s3);
    1875     1565729 :       tes = gadd(bernfrac(i), gdivgunextu(gmul(s5,tes), i+1));
    1876     1565729 :       if (gc_needed(av2,3))
    1877             :       {
    1878           0 :         if(DEBUGMEM>1) pari_warn(warnmem,"czeta i = %ld", i);
    1879           0 :         (void)gc_all(av2,3, &tes,&s5,&s4);
    1880             :       }
    1881             :     }
    1882       18445 :     u = gmul(gmul(tes,invn2), gmul2n(s2, -1));
    1883       18445 :     tes = gmulsg(nn, gaddsg(1, u));
    1884             :   }
    1885       18445 :   if (DEBUGLEVEL>2) timer_printf(&T,"Bernoulli sum");
    1886             :   /* y += tes n^(-s) / (s-1) */
    1887       18445 :   y = gadd(y, gmul(tes, gdiv(ns, gsubgs(s,1))));
    1888       18445 :   if (funeq_factor) y = gmul(y, funeq_factor);
    1889       18445 :   set_avma(av); return affc_fixlg(y,res);
    1890             : }
    1891             : /* v a t_VEC/t_COL; is v[i] = a + b (i-1) for some a,b ? */
    1892             : int
    1893          42 : RgV_is_arithprog(GEN v, GEN *a, GEN *b)
    1894             : {
    1895          42 :   pari_sp av = avma, av2;
    1896          42 :   long i, n = lg(v)-1;
    1897          42 :   if (n == 0) { *a = *b = gen_0; return 1; }
    1898          42 :   *a = gel(v,1);
    1899          42 :   if (n == 1) { * b = gen_0; return 1; }
    1900          42 :   *b = gsub(gel(v,2), *a); av2 = avma;
    1901          77 :   for (i = 2; i < n; i++)
    1902          35 :     if (!gequal(*b, gsub(gel(v,i+1), gel(v,i)))) return gc_int(av,0);
    1903          42 :   return gc_int(av2,1);
    1904             : }
    1905             : 
    1906             : GEN
    1907       23583 : gzeta(GEN x, long prec)
    1908             : {
    1909       23583 :   pari_sp av = avma;
    1910             :   GEN y;
    1911       23583 :   if (gequal1(x)) pari_err_DOMAIN("zeta", "argument", "=", gen_1, x);
    1912       23548 :   switch(typ(x))
    1913             :   {
    1914        4774 :     case t_INT:
    1915        4774 :       if (is_bigint(x))
    1916             :       {
    1917          21 :         if (signe(x) > 0) return real_1(prec);
    1918          14 :         if (mod2(x) == 0) return real_0(prec);
    1919           7 :         pari_err_OVERFLOW("zeta [large negative argument]");
    1920             :       }
    1921        4753 :       return szeta(itos(x),prec);
    1922       18515 :     case t_REAL: case t_COMPLEX: return czeta(x,prec);
    1923          14 :     case t_PADIC: return Qp_zeta(x);
    1924          49 :     case t_VEC: case t_COL:
    1925             :     {
    1926             :       GEN a, b;
    1927          49 :       long n = lg(x) - 1;
    1928          49 :       if (n > 1 && RgV_is_arithprog(x, &b, &a))
    1929             :       {
    1930          42 :         if (!is_real_t(typ(a)) || !is_real_t(typ(b))
    1931          35 :             || gcmpgs(gel(x,1), 1) <= 0
    1932          42 :             || gcmpgs(gel(x,n), 1) <= 0) { set_avma(av); break; }
    1933          21 :         a = veczeta(a, b, n, prec);
    1934          21 :         settyp(a, typ(x)); return a;
    1935             :       }
    1936             :     }
    1937             :     default:
    1938         203 :       if (!(y = toser_i(x))) break;
    1939          35 :       if (gequal1(y))
    1940           7 :         pari_err_DOMAIN("zeta", "argument", "=", gen_1, y);
    1941          28 :       return gc_upto(av, lfun(gen_1,y,prec2nbits(prec)));
    1942             :   }
    1943         189 :   return trans_eval("zeta",gzeta,x,prec);
    1944             : }
    1945             : 
    1946             : /***********************************************************************/
    1947             : /**                                                                   **/
    1948             : /**                    FONCTIONS POLYLOGARITHME                       **/
    1949             : /**                                                                   **/
    1950             : /***********************************************************************/
    1951             : 
    1952             : /* smallish k such that bernbitprec(K) > bit + Kdz, K = 2k+4 */
    1953             : static long
    1954          21 : get_k(double dz, long bit)
    1955             : {
    1956             :   long a, b;
    1957          21 :   for (b = 128;; b <<= 1)
    1958          21 :     if (bernbitprec(b) > bit + b*dz) break;
    1959          21 :   if (b == 128) return 128;
    1960           0 :   a = b >> 1;
    1961           0 :   while (b - a > 64)
    1962             :   {
    1963           0 :     long c = (a+b) >> 1;
    1964           0 :     if (bernbitprec(c) > bit + c*dz) b = c; else a = c;
    1965             :   }
    1966           0 :   return b >> 1;
    1967             : }
    1968             : 
    1969             : /* m >= 2. Validity domain |log x| < 2*Pi, contains log |x| < 5.44,
    1970             :  * Li_m(x = e^z) = sum_{n >= 0} zeta(m-n) z^n / n!
    1971             :  *    with zeta(1) := H_{m-1} - log(-z) */
    1972             : static GEN
    1973          21 : cxpolylog(long m, GEN x, long prec)
    1974             : {
    1975             :   long li, n, k, ksmall, real;
    1976             :   GEN vz, z, Z, h, q, s, S;
    1977             :   pari_sp av;
    1978             :   double dz;
    1979             :   pari_timer T;
    1980             : 
    1981          21 :   if (gequal1(x)) return szeta(m,prec);
    1982             :   /* x real <= 1 ==> Li_m(x) real */
    1983          21 :   real = (typ(x) == t_REAL && (expo(x) < 0 || signe(x) <= 0));
    1984             : 
    1985          21 :   vz = constzeta(m, prec);
    1986          21 :   z = glog(x,prec);
    1987             :   /* n = 0 */
    1988          21 :   q = gen_1; s = gel(vz, m);
    1989          28 :   for (n=1; n < m-1; n++)
    1990             :   {
    1991           7 :     q = gdivgu(gmul(q,z),n);
    1992           7 :     s = gadd(s, gmul(gel(vz,m-n), real? real_i(q): q));
    1993             :   }
    1994             :   /* n = m-1 */
    1995          21 :     q = gdivgu(gmul(q,z),n); /* multiply by "zeta(1)" */
    1996          21 :     h = gmul(q, gsub(harmonic(m-1), glog(gneg_i(z),prec)));
    1997          21 :     s = gadd(s, real? real_i(h): h);
    1998             :   /* n = m */
    1999          21 :     q = gdivgu(gmul(q,z),m);
    2000          21 :     s = gadd(s, gdivgs(real? real_i(q): q, -2)); /* zeta(0) = -1/2 */
    2001             :   /* n = m+1 */
    2002          21 :     q = gdivgu(gmul(q,z),m+1); /* = z^(m+1) / (m+1)! */
    2003          21 :     s = gadd(s, gdivgs(real? real_i(q): q, -12)); /* zeta(-1) = -1/12 */
    2004             : 
    2005          21 :   li = -(prec2nbits(prec)+1);
    2006          21 :   if (DEBUGLEVEL) timer_start(&T);
    2007          21 :   dz = dbllog2(z) - log2PI; /*  ~ log2(|z|/2Pi) */
    2008             :   /* sum_{k >= 1} zeta(-1-2k) * z^(2k+m+1) / (2k+m+1)!
    2009             :    * = 2 z^(m-1) sum_{k >= 1} zeta(2k+2) * Z^(k+1) / (2k+2)..(2k+1+m)), where
    2010             :    * Z = -(z/2Pi)^2. Stop at 2k = (li - (m-1)*Lz - m) / dz, Lz = log2 |z| */
    2011             :   /* We cut the sum in two: small values of k first */
    2012          21 :   Z = gsqr(z); av = avma;
    2013          21 :   ksmall = get_k(dz, prec2nbits(prec));
    2014          21 :   constbern(ksmall);
    2015         469 :   for(k = 1; k < ksmall; k++)
    2016             :   {
    2017         469 :     GEN t = q = gdivgunextu(gmul(q,Z), 2*k+m); /* z^(2k+m+1)/(2k+m+1)! */
    2018         469 :     if (real) t = real_i(t);
    2019         469 :     t = gmul(t, gdivgu(bernfrac(2*k+2), 2*k+2)); /* - t * zeta(1-(2k+2)) */
    2020         469 :     s = gsub(s, t);
    2021         469 :     if (gexpo(t)  < li) return s;
    2022             :     /* large values ? */
    2023         448 :     if ((k & 0x1ff) == 0) (void)gc_all(av, 2, &s, &q);
    2024             :   }
    2025           0 :   if (DEBUGLEVEL>2) timer_printf(&T, "polylog: small k <= %ld", k);
    2026           0 :   Z = gneg(gsqr(gdiv(z, Pi2n(1,prec))));
    2027           0 :   q = gmul(gpowgs(z, m-1), gpowgs(Z, k+1)); /* Z^(k+1) * z^(m-1) */
    2028           0 :   S = gen_0; av = avma;
    2029           0 :   for(;; k++)
    2030           0 :   {
    2031           0 :     GEN t = q;
    2032             :     long b;
    2033           0 :     if (real) t = real_i(t);
    2034           0 :     b = prec + gexpo(t) / BITS_IN_LONG; /* decrease accuracy */
    2035           0 :     if (b == 2) break;
    2036             :     /* t * zeta(2k+2) / (2k+2)..(2k+1+m) */
    2037           0 :     t = gdiv(t, mulri(inv_szeta_euler(2*k+2, b),
    2038           0 :                       mulu_interval(2*k+2, 2*k+1+m)));
    2039           0 :     S = gadd(S, t); if (gexpo(t)  < li) break;
    2040           0 :     q = gmul(q, Z);
    2041           0 :     if ((k & 0x1ff) == 0) (void)gc_all(av, 2, &S, &q);
    2042             :   }
    2043           0 :   if (DEBUGLEVEL>2) timer_printf(&T, "polylog: large k <= %ld", k);
    2044           0 :   return gadd(s, gmul2n(S,1));
    2045             : }
    2046             : 
    2047             : static GEN
    2048          42 : Li1(GEN x, long prec) { return gneg(glog(gsubsg(1, x), prec)); }
    2049             : 
    2050             : static GEN
    2051         203 : polylog(long m, GEN x, long prec)
    2052             : {
    2053             :   long l, e, i, G, sx;
    2054             :   pari_sp av, av1;
    2055             :   GEN X, Xn, z, p1, p2, y, res;
    2056             : 
    2057         203 :   if (m < 0) pari_err_DOMAIN("polylog", "index", "<", gen_0, stoi(m));
    2058         203 :   if (!m) return mkfrac(gen_m1,gen_2);
    2059         203 :   if (gequal0(x)) return gcopy(x);
    2060         203 :   if (m==1) { av = avma; return gc_upto(av, Li1(x, prec)); }
    2061             : 
    2062         168 :   l = precision(x);
    2063         168 :   if (!l) l = prec; else prec = l;
    2064         168 :   res = cgetc(l); av = avma;
    2065         168 :   x = gtofp(x, l+EXTRAPREC64);
    2066         168 :   e = gexpo(gnorm(x));
    2067         168 :   if (!e || e == -1) {
    2068          21 :     y = cxpolylog(m,x,prec);
    2069          21 :     set_avma(av); return affc_fixlg(y, res);
    2070             :   }
    2071         147 :   X = (e > 0)? ginv(x): x;
    2072         147 :   G = -prec2nbits(l);
    2073         147 :   av1 = avma;
    2074         147 :   y = Xn = X;
    2075       68159 :   for (i=2; ; i++)
    2076             :   {
    2077       68159 :     Xn = gmul(X,Xn); p2 = gdiv(Xn,powuu(i,m));
    2078       68159 :     y = gadd(y,p2);
    2079       68159 :     if (gexpo(p2) <= G) break;
    2080             : 
    2081       68012 :     if (gc_needed(av1,1))
    2082             :     {
    2083           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"polylog");
    2084           0 :       (void)gc_all(av1,2, &y, &Xn);
    2085             :     }
    2086             :   }
    2087         147 :   if (e < 0) { set_avma(av); return affc_fixlg(y, res); }
    2088             : 
    2089          28 :   sx = gsigne(imag_i(x));
    2090          28 :   if (!sx)
    2091             :   {
    2092          28 :     if (m&1) sx = gsigne(gsub(gen_1, real_i(x)));
    2093          21 :     else     sx = - gsigne(real_i(x));
    2094             :   }
    2095          28 :   z = divri(mppi(l), mpfact(m-1)); setsigne(z, sx);
    2096          28 :   z = mkcomplex(gen_0, z);
    2097             : 
    2098          28 :   if (m == 2)
    2099             :   { /* same formula as below, written more efficiently */
    2100          21 :     y = gneg_i(y);
    2101          21 :     if (typ(x) == t_REAL && signe(x) < 0)
    2102           7 :       p1 = logr_abs(x);
    2103             :     else
    2104          14 :       p1 = gsub(glog(x,l), z);
    2105          21 :     p1 = gmul2n(gsqr(p1), -1); /* = (log(-x))^2 / 2 */
    2106             : 
    2107          21 :     p1 = gadd(p1, divru(sqrr(mppi(l)), 6));
    2108          21 :     p1 = gneg_i(p1);
    2109             :   }
    2110             :   else
    2111             :   {
    2112           7 :     GEN logx = glog(x,l), logx2 = gsqr(logx), vz = constzeta(m, l);
    2113           7 :     p1 = mkfrac(gen_m1,gen_2);
    2114          14 :     for (i = m-2; i >= 0; i -= 2)
    2115           7 :       p1 = gadd(gel(vz, m-i), gmul(p1, gdivgunextu(logx2, i+1)));
    2116           7 :     if (m&1) p1 = gmul(logx,p1); else y = gneg_i(y);
    2117           7 :     p1 = gadd(gmul2n(p1,1), gmul(z,gpowgs(logx,m-1)));
    2118           7 :     if (typ(x) == t_REAL && signe(x) < 0) p1 = real_i(p1);
    2119             :   }
    2120          28 :   y = gadd(y,p1);
    2121          28 :   set_avma(av); return affc_fixlg(y, res);
    2122             : }
    2123             : static GEN
    2124         119 : RIpolylog(long m, GEN x, long real, long prec)
    2125             : {
    2126         119 :   GEN y = polylog(m, x, prec);
    2127         119 :   return real? real_i(y): imag_i(y);
    2128             : }
    2129             : GEN
    2130          21 : dilog(GEN x, long prec) { return gpolylog(2, x, prec); }
    2131             : 
    2132             : /* x a floating point number, t_REAL or t_COMPLEX of t_REAL */
    2133             : static GEN
    2134          42 : logabs(GEN x)
    2135             : {
    2136             :   GEN y;
    2137          42 :   if (typ(x) == t_COMPLEX)
    2138             :   {
    2139           7 :     y = logr_abs( cxnorm(x) );
    2140           7 :     shiftr_inplace(y, -1);
    2141             :   } else
    2142          35 :     y = logr_abs(x);
    2143          42 :   return y;
    2144             : }
    2145             : 
    2146             : static GEN
    2147          21 : polylogD(long m, GEN x, long flag, long prec)
    2148             : {
    2149          21 :   long fl = 0, k, l, m2;
    2150             :   pari_sp av;
    2151             :   GEN p1, p2, y;
    2152             : 
    2153          21 :   if (gequal0(x)) return gcopy(x);
    2154          21 :   m2 = m&1;
    2155          21 :   if (gequal1(x) && m>=2) return m2? szeta(m,prec): gen_0;
    2156          21 :   av = avma; l = precision(x);
    2157          21 :   if (!l) { l = prec; x = gtofp(x,l); }
    2158          21 :   p1 = logabs(x);
    2159          21 :   if (signe(p1) > 0) { x = ginv(x); fl = !m2; } else setabssign(p1);
    2160             :   /* |x| <= 1, p1 = - log|x| >= 0 */
    2161          21 :   p2 = gen_1;
    2162          21 :   y = RIpolylog(m, x, m2, l);
    2163          84 :   for (k = 1; k < m; k++)
    2164             :   {
    2165          63 :     GEN t = RIpolylog(m-k, x, m2, l);
    2166          63 :     p2 = gdivgu(gmul(p2,p1), k); /* (-log|x|)^k / k! */
    2167          63 :     y = gadd(y, gmul(p2, t));
    2168             :   }
    2169          21 :   if (m2)
    2170             :   {
    2171          14 :     p1 = flag? gdivgs(p1, -2*m): gdivgs(logabs(gsubsg(1,x)), m);
    2172          14 :     y = gadd(y, gmul(p2, p1));
    2173             :   }
    2174          21 :   if (fl) y = gneg(y);
    2175          21 :   return gc_upto(av, y);
    2176             : }
    2177             : 
    2178             : static GEN
    2179          14 : polylogP(long m, GEN x, long prec)
    2180             : {
    2181          14 :   long fl = 0, k, l, m2;
    2182             :   pari_sp av;
    2183             :   GEN p1,y;
    2184             : 
    2185          14 :   if (gequal0(x)) return gcopy(x);
    2186          14 :   m2 = m&1;
    2187          14 :   if (gequal1(x) && m>=2) return m2? szeta(m,prec): gen_0;
    2188          14 :   av = avma; l = precision(x);
    2189          14 :   if (!l) { l = prec; x = gtofp(x,l); }
    2190          14 :   p1 = logabs(x);
    2191          14 :   if (signe(p1) > 0) { x = ginv(x); fl = !m2; setsigne(p1, -1); }
    2192             :   /* |x| <= 1 */
    2193          14 :   y = RIpolylog(m, x, m2, l);
    2194          14 :   if (m==1)
    2195             :   {
    2196           7 :     shiftr_inplace(p1, -1); /* log |x| / 2 */
    2197           7 :     y = gadd(y, p1);
    2198             :   }
    2199             :   else
    2200             :   { /* m >= 2, \sum_{0 <= k <= m} 2^k B_k/k! (log |x|)^k Li_{m-k}(x),
    2201             :        with Li_0(x) := -1/2 */
    2202           7 :     GEN u, t = RIpolylog(m-1, x, m2, l);
    2203           7 :     u = gneg_i(p1); /* u = 2 B1 log |x| */
    2204           7 :     y = gadd(y, gmul(u, t));
    2205           7 :     if (m > 2)
    2206             :     {
    2207             :       GEN p2;
    2208           7 :       shiftr_inplace(p1, 1); /* 2log|x| <= 0 */
    2209           7 :       constbern(m>>1);
    2210           7 :       p1 = sqrr(p1);
    2211           7 :       p2 = shiftr(p1,-1);
    2212          21 :       for (k = 2; k < m; k += 2)
    2213             :       {
    2214          14 :         if (k > 2) p2 = gdivgunextu(gmul(p2,p1),k-1); /* 2^k/k! log^k |x|*/
    2215          14 :         t = RIpolylog(m-k, x, m2, l);
    2216          14 :         u = gmul(p2, bernfrac(k));
    2217          14 :         y = gadd(y, gmul(u, t));
    2218             :       }
    2219             :     }
    2220             :   }
    2221          14 :   if (fl) y = gneg(y);
    2222          14 :   return gc_upto(av, y);
    2223             : }
    2224             : 
    2225             : static GEN
    2226         175 : gpolylog_i(void *E, GEN x, long prec)
    2227             : {
    2228         175 :   pari_sp av = avma;
    2229         175 :   long i, n, v, m = (long)E;
    2230             :   GEN a, y;
    2231             : 
    2232         175 :   if (m <= 0)
    2233             :   {
    2234          28 :     a = gmul(x, poleval(eulerianpol(-m, 0), x));
    2235          28 :     return gc_upto(av, gdiv(a, gpowgs(gsubsg(1, x), 1-m)));
    2236             :   }
    2237         147 :   switch(typ(x))
    2238             :   {
    2239          84 :     case t_REAL: case t_COMPLEX: return polylog(m,x,prec);
    2240           7 :     case t_INTMOD: case t_PADIC: pari_err_IMPL( "padic polylogarithm");
    2241          56 :     default:
    2242          56 :       av = avma; if (!(y = toser_i(x))) break;
    2243          21 :       if (!m) { set_avma(av); return mkfrac(gen_m1,gen_2); }
    2244          21 :       if (m==1) return gc_upto(av, Li1(y, prec));
    2245          21 :       if (gequal0(y)) return gc_GEN(av, y);
    2246          21 :       v = valser(y);
    2247          21 :       if (v < 0) pari_err_DOMAIN("polylog","valuation", "<", gen_0, x);
    2248          14 :       if (v > 0) {
    2249           7 :         n = (lg(y)-3 + v) / v;
    2250           7 :         a = zeroser(varn(y), lg(y)-2);
    2251          35 :         for (i=n; i>=1; i--)
    2252          28 :           a = gmul(y, gadd(a, powis(utoipos(i),-m)));
    2253             :       } else { /* v == 0 */
    2254           7 :         long vy = varn(y);
    2255           7 :         GEN a0 = polcoef_i(y, 0, -1), t = gdiv(derivser(y), y);
    2256           7 :         a = Li1(y, prec);
    2257          14 :         for (i=2; i<=m; i++)
    2258           7 :           a = gadd(gpolylog(i, a0, prec), integ(gmul(t, a), vy));
    2259             :       }
    2260          14 :       return gc_upto(av, a);
    2261             :   }
    2262          35 :   return trans_evalgen("polylog", E, gpolylog_i, x, prec);
    2263             : }
    2264             : GEN
    2265         133 : gpolylog(long m, GEN x, long prec) { return gpolylog_i((void*)m, x, prec); }
    2266             : 
    2267             : GEN
    2268         147 : polylog0(long m, GEN x, long flag, long prec)
    2269             : {
    2270         147 :   switch(flag)
    2271             :   {
    2272         105 :     case 0: return gpolylog(m,x,prec);
    2273          14 :     case 1: return polylogD(m,x,0,prec);
    2274           7 :     case 2: return polylogD(m,x,1,prec);
    2275          14 :     case 3: return polylogP(m,x,prec);
    2276           7 :     default: pari_err_FLAG("polylog");
    2277             :   }
    2278             :   return NULL; /* LCOV_EXCL_LINE */
    2279             : }

Generated by: LCOV version 1.16