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.16.2 lcov report (development 29115-f22e516b23) Lines: 1739 1866 93.2 %
Date: 2024-03-19 08:06:49 Functions: 126 129 97.7 %
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 = gerepileupto(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 gerepileupto(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 gerepileupto(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 gerepileupto(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 gerepileupto(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 gerepileupto(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 gerepileupto(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 gerepileupto(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 gerepileupto(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)) gerepileall(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)) gerepileall(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 :       gerepileall(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 = gerepileuptoleaf(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 = gerepileupto(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 gerepileupto(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 gerepilecopy(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 gerepilecopy(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 gerepilecopy(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 gerepileupto(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 gerepileupto(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 gerepileupto(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 gerepilecopy(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       16639 : gamma_use_asymp(GEN x, long b)
     720             : {
     721             :   long e;
     722       16639 :   if (is_real_t(typ(x)))
     723             :   {
     724       12607 :     pari_sp av = avma;
     725       12607 :     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)) gerepileall(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 gerepileuptoleaf(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)) gerepileall(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 gerepileupto(av, gmul(S, ix));
     782             : }
     783             : 
     784             : /* eint1(x) = incgam(0, x); typ(x) = t_REAL, x > 0 */
     785             : static GEN
     786        6398 : eint1p(GEN x, GEN expx)
     787             : {
     788             :   pari_sp av;
     789        6398 :   long l = realprec(x), bit = prec2nbits(l), prec, i;
     790             :   double mx;
     791             :   GEN z, S, t, H, run;
     792             : 
     793        6398 :   if (gamma_use_asymp(x, bit)
     794          28 :       && (z = eint1r_asymp(x, expx, l))) return z;
     795        6370 :   mx = rtodbl(x);
     796        6370 :   prec = l + nbits2extraprec((mx+log(mx))/M_LN2 + 10);
     797        6370 :   bit = prec2nbits(prec);
     798        6370 :   run = real_1(prec); x = rtor(x, prec);
     799        6370 :   av = avma; S = z = t = H = run;
     800      608487 :   for (i = 2; expo(S) - expo(t) <= bit; i++)
     801             :   {
     802      602117 :     H = addrr(H, divru(run,i)); /* H = sum_{k<=i} 1/k */
     803      602117 :     z = divru(mulrr(x,z), i);   /* z = x^(i-1)/i! */
     804      602117 :     t = mulrr(z, H); S = addrr(S, t);
     805      602117 :     if ((i & 0x1ff) == 0) gerepileall(av, 4, &z,&t,&S,&H);
     806             :   }
     807        6370 :   return subrr(mulrr(x, divrr(S,expx? expx: mpexp(x))),
     808             :                addrr(mplog(x), mpeuler(prec)));
     809             : }
     810             : /* eint1(x) = incgam(0, x); typ(x) = t_REAL, x < 0
     811             :  * rewritten from code contributed by Manfred Radimersky */
     812             : static GEN
     813         140 : eint1m(GEN x, GEN expx)
     814             : {
     815         140 :   GEN p1, q, S, y, z = cgetg(3, t_COMPLEX);
     816         140 :   long l  = realprec(x), n  = prec2nbits(l), j;
     817         140 :   pari_sp av = avma;
     818             : 
     819         140 :   y  = rtor(x, l + EXTRAPREC64); setsigne(y,1); /* |x| */
     820         140 :   if (gamma_use_asymp(y, n))
     821             :   { /* ~eint1_asymp: asymptotic expansion */
     822          14 :     p1 = q = invr(y); S = addrs(q, 1);
     823         560 :     for (j = 2; expo(q) >= -n; j++) {
     824         546 :       q = mulrr(q, mulru(p1, j));
     825         546 :       S = addrr(S, q);
     826             :     }
     827          14 :     y  = mulrr(p1, expx? divrr(S, expx): mulrr(S, mpexp(y)));
     828             :   }
     829             :   else
     830             :   {
     831         126 :     p1 = q = S = y;
     832       24248 :     for (j = 2; expo(q) - expo(S) >= -n; j++) {
     833       24122 :       p1 = mulrr(y, divru(p1, j)); /* (-x)^j/j! */
     834       24122 :       q = divru(p1, j);
     835       24122 :       S = addrr(S, q);
     836             :     }
     837         126 :     y  = addrr(S, addrr(logr_abs(x), mpeuler(l)));
     838             :   }
     839         140 :   y = gerepileuptoleaf(av, y); togglesign(y);
     840         140 :   gel(z, 1) = y;
     841         140 :   y = mppi(l); setsigne(y, -1);
     842         140 :   gel(z, 2) = y; return z;
     843             : }
     844             : 
     845             : /* real(z*log(z)-z), z = x+iy */
     846             : static double
     847        8372 : mygamma(double x, double y)
     848             : {
     849        8372 :   if (x == 0.) return -(M_PI/2)*fabs(y);
     850        8372 :   return (x/2)*log(x*x+y*y)-x-y*atan(y/x);
     851             : }
     852             : 
     853             : /* x^s exp(-x) */
     854             : static GEN
     855       10843 : expmx_xs(GEN s, GEN x, GEN logx, long prec)
     856             : {
     857             :   GEN z;
     858       10843 :   long ts = typ(s);
     859       10843 :   if (ts == t_INT || (ts == t_FRAC && absequaliu(gel(s,2), 2)))
     860        5264 :     z = gmul(gexp(gneg(x), prec), gpow(x, s, prec));
     861             :   else
     862        5579 :     z = gexp(gsub(gmul(s, logx? logx: glog(x,prec+EXTRAPREC64)), x), prec);
     863       10843 :   return z;
     864             : }
     865             : 
     866             : /* Not yet: doesn't work at low accuracy
     867             : #define INCGAM_CF
     868             : */
     869             : 
     870             : #ifdef INCGAM_CF
     871             : /* Is s very close to a nonpositive integer ? */
     872             : static int
     873             : isgammapole(GEN s, long bitprec)
     874             : {
     875             :   pari_sp av = avma;
     876             :   GEN t = imag_i(s);
     877             :   long e, b = bitprec - 10;
     878             : 
     879             :   if (gexpo(t) > - b) return 0;
     880             :   s = real_i(s);
     881             :   if (gsigne(s) > 0 && gexpo(s) > -b) return 0;
     882             :   (void)grndtoi(s, &e); return gc_bool(av, e < -b);
     883             : }
     884             : 
     885             : /* incgam using the continued fraction. x a t_REAL or t_COMPLEX, mx ~ |x|.
     886             :  * Assume precision(s), precision(x) >= prec */
     887             : static GEN
     888             : incgam_cf(GEN s, GEN x, double mx, long prec)
     889             : {
     890             :   GEN ms, y, S;
     891             :   long n, i, j, LS, bitprec = prec2nbits(prec);
     892             :   double rs, is, m;
     893             : 
     894             :   if (typ(s) == t_COMPLEX)
     895             :   {
     896             :     rs = gtodouble(gel(s,1));
     897             :     is = gtodouble(gel(s,2));
     898             :   }
     899             :   else
     900             :   {
     901             :     rs = gtodouble(s);
     902             :     is = 0.;
     903             :   }
     904             :   if (isgammapole(s, bitprec)) LS = 0;
     905             :   else
     906             :   {
     907             :     double bit,  LGS = mygamma(rs,is);
     908             :     LS = LGS <= 0 ? 0: ceil(LGS);
     909             :     bit = (LGS - (rs-1)*log(mx) + mx)/M_LN2;
     910             :     if (bit > 0)
     911             :     {
     912             :       prec += nbits2extraprec((long)bit);
     913             :       x = gtofp(x, prec);
     914             :       if (isinexactreal(s)) s = gtofp(s, prec);
     915             :     }
     916             :   }
     917             :   /* |ln(2*gamma(s)*sin(s*Pi))| <= ln(2) + |lngamma(s)| + |Im(s)*Pi|*/
     918             :   m = bitprec*M_LN2 + LS + M_LN2 + fabs(is)*M_PI + mx;
     919             :   if (rs < 1) m += (1 - rs)*log(mx);
     920             :   m /= 4;
     921             :   n = (long)(1 + m*m/mx);
     922             :   y = expmx_xs(gsubgs(s,1), x, NULL, prec);
     923             :   if (rs >= 0 && bitprec >= 512)
     924             :   {
     925             :     GEN A = cgetg(n+1, t_VEC), B = cgetg(n+1, t_VEC);
     926             :     ms = gsubsg(1, s);
     927             :     for (j = 1; j <= n; ++j)
     928             :     {
     929             :       gel(A,j) = ms;
     930             :       gel(B,j) = gmulsg(j, gsubgs(s,j));
     931             :       ms = gaddgs(ms, 2);
     932             :     }
     933             :     S = contfraceval_inv(mkvec2(A,B), x, -1);
     934             :   }
     935             :   else
     936             :   {
     937             :     GEN x_s = gsub(x, s);
     938             :     pari_sp av2 = avma;
     939             :     S = gdiv(gsubgs(s,n), gaddgs(x_s,n<<1));
     940             :     for (i=n-1; i >= 1; i--)
     941             :     {
     942             :       S = gdiv(gsubgs(s,i), gadd(gaddgs(x_s,i<<1),gmulsg(i,S)));
     943             :       if (gc_needed(av2,3))
     944             :       {
     945             :         if(DEBUGMEM>1) pari_warn(warnmem,"incgam_cf");
     946             :         S = gerepileupto(av2, S);
     947             :       }
     948             :     }
     949             :     S = gaddgs(S,1);
     950             :   }
     951             :   return gmul(y, S);
     952             : }
     953             : #endif
     954             : 
     955             : static double
     956        6419 : findextraincgam(GEN s, GEN x)
     957             : {
     958        6419 :   double sig = gtodouble(real_i(s)), t = gtodouble(imag_i(s));
     959        6419 :   double xr = gtodouble(real_i(x)), xi = gtodouble(imag_i(x));
     960        6419 :   double exd = 0., Nx = xr*xr + xi*xi, D = Nx - t*t;
     961             :   long n;
     962             : 
     963        6419 :   if (xr < 0)
     964             :   {
     965         833 :     long ex = gexpo(x);
     966         833 :     if (ex > 0 && ex > gexpo(s)) exd = sqrt(Nx)*log(Nx)/2; /* |x| log |x| */
     967             :   }
     968        6419 :   if (D <= 0.) return exd;
     969        4977 :   n = (long)(sqrt(D)-sig);
     970        4977 :   if (n <= 0) return exd;
     971        1841 :   return maxdd(exd, (n*log(Nx)/2 - mygamma(sig+n, t) + mygamma(sig, t)) / M_LN2);
     972             : }
     973             : 
     974             : /* use exp(-x) * (x^s/s) * sum_{k >= 0} x^k / prod(i=1, k, s+i) */
     975             : static GEN
     976        6426 : incgamc_i(GEN s, GEN x, long *ptexd, long prec)
     977             : {
     978             :   GEN S, t, y;
     979             :   long l, n, i, exd;
     980        6426 :   pari_sp av = avma, av2;
     981             : 
     982        6426 :   if (gequal0(x))
     983             :   {
     984           7 :     if (ptexd) *ptexd = 0.;
     985           7 :     return gtofp(x, prec);
     986             :   }
     987        6419 :   l = precision(x);
     988        6419 :   if (!l) l = prec;
     989        6419 :   n = -prec2nbits(l)-1;
     990        6419 :   exd = (long)findextraincgam(s, x);
     991        6419 :   if (ptexd) *ptexd = exd;
     992        6419 :   if (exd > 0)
     993             :   {
     994        1666 :     long p = l + nbits2extraprec(exd);
     995        1666 :     x = gtofp(x, p);
     996        1666 :     if (isinexactreal(s)) s = gtofp(s, p);
     997             :   }
     998        4753 :   else x = gtofp(x, l+EXTRAPREC64);
     999        6419 :   av2 = avma;
    1000        6419 :   S = gdiv(x, gaddsg(1,s));
    1001        6419 :   t = gaddsg(1, S);
    1002      770875 :   for (i=2; gexpo(S) >= n; i++)
    1003             :   {
    1004      764456 :     S = gdiv(gmul(x,S), gaddsg(i,s)); /* x^i / ((s+1)...(s+i)) */
    1005      764456 :     t = gadd(S,t);
    1006      764456 :     if (gc_needed(av2,3))
    1007             :     {
    1008           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"incgamc");
    1009           0 :       gerepileall(av2, 2, &S, &t);
    1010             :     }
    1011             :   }
    1012        6419 :   y = expmx_xs(s, x, NULL, prec);
    1013        6419 :   return gerepileupto(av, gmul(gdiv(y,s), t));
    1014             : }
    1015             : 
    1016             : GEN
    1017        2226 : incgamc(GEN s, GEN x, long prec)
    1018        2226 : { return incgamc_i(s, x, NULL, prec); }
    1019             : 
    1020             : /* incgamma using asymptotic expansion:
    1021             :  *   exp(-x)x^(s-1)(1 + (s-1)/x + (s-1)(s-2)/x^2 + ...) */
    1022             : static GEN
    1023        2716 : incgam_asymp(GEN s, GEN x, long prec)
    1024             : {
    1025        2716 :   pari_sp av = avma, av2;
    1026             :   GEN S, q, cox, invx;
    1027        2716 :   long oldeq = LONG_MAX, eq, esx, j;
    1028        2716 :   int flint = (typ(s) == t_INT && signe(s) > 0);
    1029             : 
    1030        2716 :   x = gtofp(x,prec+EXTRAPREC64);
    1031        2716 :   invx = ginv(x);
    1032        2716 :   esx = -prec2nbits(prec);
    1033        2716 :   av2 = avma;
    1034        2716 :   q = gmul(gsubgs(s,1), invx);
    1035        2716 :   S = gaddgs(q, 1);
    1036        2716 :   for (j = 2;; j++)
    1037             :   {
    1038      123879 :     eq = gexpo(q); if (eq < esx) break;
    1039      121282 :     if (!flint && (j & 3) == 0)
    1040             :     { /* guard against divergence */
    1041       15778 :       if (eq > oldeq) return gc_NULL(av); /* regressing, abort */
    1042       15659 :       oldeq = eq;
    1043             :     }
    1044      121163 :     q = gmul(q, gmul(gsubgs(s,j), invx));
    1045      121163 :     S = gadd(S, q);
    1046      121163 :     if (gc_needed(av2, 1)) gerepileall(av2, 2, &S, &q);
    1047             :   }
    1048        2597 :   if (DEBUGLEVEL > 2) err_printf("incgam: using asymp\n");
    1049        2597 :   cox = expmx_xs(gsubgs(s,1), x, NULL, prec);
    1050        2597 :   return gerepileupto(av, gmul(cox, S));
    1051             : }
    1052             : 
    1053             : /* gasx = incgam(s-n,x). Compute incgam(s,x)
    1054             :  * = (s-1)(s-2)...(s-n)gasx + exp(-x)x^(s-1) *
    1055             :  *   (1 + (s-1)/x + ... + (s-1)(s-2)...(s-n+1)/x^(n-1)) */
    1056             : static GEN
    1057         546 : incgam_asymp_partial(GEN s, GEN x, GEN gasx, long n, long prec)
    1058             : {
    1059             :   pari_sp av;
    1060         546 :   GEN S, q, cox, invx, s1 = gsubgs(s, 1), sprod;
    1061             :   long j;
    1062         546 :   cox = expmx_xs(s1, x, NULL, prec);
    1063         546 :   if (n == 1) return gadd(cox, gmul(s1, gasx));
    1064         546 :   invx = ginv(x);
    1065         546 :   av = avma;
    1066         546 :   q = gmul(s1, invx);
    1067         546 :   S = gaddgs(q, 1);
    1068       52164 :   for (j = 2; j < n; j++)
    1069             :   {
    1070       51618 :     q = gmul(q, gmul(gsubgs(s, j), invx));
    1071       51618 :     S = gadd(S, q);
    1072       51618 :     if (gc_needed(av, 2)) gerepileall(av, 2, &S, &q);
    1073             :   }
    1074         546 :   sprod = gmul(gmul(q, gpowgs(x, n-1)), gsubgs(s, n));
    1075         546 :   return gadd(gmul(cox, S), gmul(sprod, gasx));
    1076             : }
    1077             : 
    1078             : /* Assume s != 0; called when Re(s) <= 1/2 */
    1079             : static GEN
    1080        2401 : incgamspec(GEN s, GEN x, GEN g, long prec)
    1081             : {
    1082        2401 :   GEN q, S, cox = gen_0, P, sk, S1, S2, S3, F3, logx, mx;
    1083        2401 :   long n, esk, E, k = itos(ground(gneg(real_i(s)))); /* >= 0 */
    1084             : 
    1085        2401 :   if (k && gexpo(x) > 0)
    1086             :   {
    1087         245 :     GEN xk = gdivgu(x, k);
    1088         245 :     long bitprec = prec2nbits(prec);
    1089         245 :     double d = (gexpo(xk) > bitprec)? bitprec*M_LN2: log(dblmodulus(xk));
    1090         245 :     d = k * (d + 1) / M_LN2;
    1091         245 :     if (d > 0) prec += nbits2extraprec((long)d);
    1092         245 :     if (isinexactreal(s)) s = gtofp(s, prec);
    1093             :   }
    1094        2401 :   x = gtofp(x, maxss(precision(x), prec) + EXTRAPREC64);
    1095        2401 :   sk = gaddgs(s, k); /* |Re(sk)| <= 1/2 */
    1096        2401 :   logx = glog(x, prec);
    1097        2401 :   mx = gneg(x);
    1098        2401 :   if (k == 0) { S = gen_0; P = gen_1; }
    1099             :   else
    1100             :   {
    1101             :     long j;
    1102         854 :     q = ginv(s); S = q; P = s;
    1103       16926 :     for (j = 1; j < k; j++)
    1104             :     {
    1105       16072 :       GEN sj = gaddgs(s, j);
    1106       16072 :       q = gmul(q, gdiv(x, sj));
    1107       16072 :       S = gadd(S, q);
    1108       16072 :       P = gmul(P, sj);
    1109             :     }
    1110         854 :     cox = expmx_xs(s, x, logx, prec); /* x^s exp(-x) */
    1111         854 :     S = gmul(S, gneg(cox));
    1112             :   }
    1113        2401 :   if (k && gequal0(sk))
    1114         175 :     return gadd(S, gdiv(eint1(x, prec), P));
    1115        2226 :   esk = gexpo(sk);
    1116        2226 :   if (esk > -7)
    1117             :   {
    1118        1015 :     GEN a, b, PG = gmul(sk, P);
    1119        1015 :     if (g) g = gmul(g, PG);
    1120        1015 :     a = incgam0(gaddgs(sk,1), x, g, prec);
    1121        1015 :     if (k == 0) cox = expmx_xs(s, x, logx, prec);
    1122        1015 :     b = gmul(gpowgs(x, k), cox);
    1123        1015 :     return gadd(S, gdiv(gsub(a, b), PG));
    1124             :   }
    1125        1211 :   E = prec2nbits(prec) + 1;
    1126        1211 :   if (gexpo(x) > 0)
    1127             :   {
    1128         420 :     long X = (long)(dblmodulus(x)/M_LN2);
    1129         420 :     prec += 2*nbits2extraprec(X);
    1130         420 :     x = gtofp(x, prec); mx = gneg(x);
    1131         420 :     logx = glog(x, prec); sk = gtofp(sk, prec);
    1132         420 :     E += X;
    1133             :   }
    1134        1211 :   if (isinexactreal(sk)) sk = gtofp(sk, prec+EXTRAPREC64);
    1135             :   /* |sk| < 2^-7 is small, guard against cancellation */
    1136        1211 :   F3 = gexpm1(gmul(sk, logx), prec);
    1137             :   /* ( gamma(1+sk) - exp(sk log(x))) ) / sk */
    1138        1211 :   S1 = gdiv(gsub(ggamma1m1(sk, prec+EXTRAPREC64), F3), sk);
    1139        1211 :   q = x; S3 = gdiv(x, gaddsg(1,sk));
    1140      255523 :   for (n = 2; gexpo(q) - gexpo(S3) > -E; n++)
    1141             :   {
    1142      254312 :     q = gmul(q, gdivgu(mx, n));
    1143      254312 :     S3 = gadd(S3, gdiv(q, gaddsg(n, sk)));
    1144             :   }
    1145        1211 :   S2 = gadd(gadd(S1, S3), gmul(F3, S3));
    1146        1211 :   return gadd(S, gdiv(S2, P));
    1147             : }
    1148             : 
    1149             : /* return |x| */
    1150             : double
    1151    11605759 : dblmodulus(GEN x)
    1152             : {
    1153    11605759 :   if (typ(x) == t_COMPLEX)
    1154             :   {
    1155     1803571 :     double a = gtodouble(gel(x,1));
    1156     1803571 :     double b = gtodouble(gel(x,2));
    1157     1803571 :     return sqrt(a*a + b*b);
    1158             :   }
    1159             :   else
    1160     9802188 :     return fabs(gtodouble(x));
    1161             : }
    1162             : 
    1163             : /* Driver routine. If g != NULL, assume that g=gamma(s,prec). */
    1164             : GEN
    1165       11550 : incgam0(GEN s, GEN x, GEN g, long prec)
    1166             : {
    1167             :   pari_sp av;
    1168             :   long E, l;
    1169             :   GEN z, rs, is;
    1170             : 
    1171       11550 :   if (gequal0(x)) return g? gcopy(g): ggamma(s,prec);
    1172       11550 :   if (gequal0(s)) return eint1(x, prec);
    1173        9744 :   l = precision(s); if (!l) l = prec;
    1174        9744 :   E = prec2nbits(l);
    1175        9744 :   if (gamma_use_asymp(x, E) ||
    1176        8323 :       (typ(s) == t_INT && signe(s) > 0 && gexpo(x) >= expi(s)))
    1177        2716 :     if ((z = incgam_asymp(s, x, l))) return z;
    1178        7147 :   av = avma; E++;
    1179        7147 :   rs = real_i(s);
    1180        7147 :   is = imag_i(s);
    1181             : #ifdef INCGAM_CF
    1182             :   /* Can one use continued fraction ? */
    1183             :   if (gequal0(is) && gequal0(imag_i(x)) && gsigne(x) > 0)
    1184             :   {
    1185             :     double sd = gtodouble(rs), LB, UB;
    1186             :     double xd = gtodouble(real_i(x));
    1187             :     if (sd > 0) {
    1188             :       LB = 15 + 0.1205*E;
    1189             :       UB = 5 + 0.752*E;
    1190             :     } else {
    1191             :       LB = -6 + 0.1205*E;
    1192             :       UB = 5 + 0.752*E + fabs(sd)/54.;
    1193             :     }
    1194             :     if (xd >= LB && xd <= UB)
    1195             :     {
    1196             :       if (DEBUGLEVEL > 2) err_printf("incgam: using continued fraction\n");
    1197             :       return gerepileupto(av, incgam_cf(s, x, xd, prec));
    1198             :     }
    1199             :   }
    1200             : #endif
    1201        7147 :   if (gsigne(rs) > 0 && gexpo(rs) >= -1)
    1202             :   { /* use complementary incomplete gamma */
    1203        4746 :     long n, egs, exd, precg, es = gexpo(s);
    1204        4746 :     if (es < 0) {
    1205         602 :       l += nbits2extraprec(-es) + 1;
    1206         602 :       x = gtofp(x, l);
    1207         602 :       if (isinexactreal(s)) s = gtofp(s, l);
    1208             :     }
    1209        4746 :     n = itos(gceil(rs));
    1210        4746 :     if (n > 100)
    1211             :     {
    1212             :       GEN gasx;
    1213         546 :       n -= 100;
    1214         546 :       if (es > 0)
    1215             :       {
    1216         546 :         es = mygamma(gtodouble(rs) - n, gtodouble(is)) / M_LN2;
    1217         546 :         if (es > 0)
    1218             :         {
    1219         546 :           l += nbits2extraprec(es);
    1220         546 :           x = gtofp(x, l);
    1221         546 :           if (isinexactreal(s)) s = gtofp(s, l);
    1222             :         }
    1223             :       }
    1224         546 :       gasx = incgam0(gsubgs(s, n), x, NULL, prec);
    1225         546 :       return gerepileupto(av, incgam_asymp_partial(s, x, gasx, n, prec));
    1226             :     }
    1227        4200 :     if (DEBUGLEVEL > 2) err_printf("incgam: using power series 1\n");
    1228             :     /* egs ~ expo(gamma(s)) */
    1229        4200 :     precg = g? precision(g): 0;
    1230        4200 :     egs = g? gexpo(g): (long)(mygamma(gtodouble(rs), gtodouble(is)) / M_LN2);
    1231        4200 :     if (egs > 0) {
    1232        1946 :       l += nbits2extraprec(egs) + 1;
    1233        1946 :       x = gtofp(x, l);
    1234        1946 :       if (isinexactreal(s)) s = gtofp(s, l);
    1235        1946 :       if (precg < l) g = NULL;
    1236             :     }
    1237        4200 :     z = incgamc_i(s, x, &exd, l);
    1238        4200 :     if (exd > 0)
    1239             :     {
    1240         896 :       l += nbits2extraprec(exd);
    1241         896 :       if (isinexactreal(s)) s = gtofp(s, l);
    1242         896 :       if (precg < l) g = NULL;
    1243             :     }
    1244             :     else
    1245             :     { /* gamma(s) negligible ? Compute to lower accuracy */
    1246        3304 :       long e = gexpo(z) - egs;
    1247        3304 :       if (e > 3)
    1248             :       {
    1249         420 :         E -= e;
    1250         420 :         if (E <= 0) g = gen_0; else if (!g) g = ggamma(s, nbits2prec(E));
    1251             :       }
    1252             :     }
    1253             :     /* worry about possible cancellation */
    1254        4200 :     if (!g) g = ggamma(s, maxss(l,precision(z)));
    1255        4200 :     return gerepileupto(av, gsub(g,z));
    1256             :   }
    1257        2401 :   if (DEBUGLEVEL > 2) err_printf("incgam: using power series 2\n");
    1258        2401 :   return gerepileupto(av, incgamspec(s, x, g, l));
    1259             : }
    1260             : 
    1261             : GEN
    1262        1106 : incgam(GEN s, GEN x, long prec) { return incgam0(s, x, NULL, prec); }
    1263             : 
    1264             : /* x a t_REAL */
    1265             : GEN
    1266        2926 : mpeint1(GEN x, GEN expx)
    1267             : {
    1268        2926 :   long s = signe(x);
    1269             :   pari_sp av;
    1270             :   GEN z;
    1271        2926 :   if (!s) pari_err_DOMAIN("eint1", "x","=",gen_0, x);
    1272        2919 :   if (s < 0) return eint1m(x, expx);
    1273        2779 :   z = cgetr(realprec(x));
    1274        2779 :   av = avma; affrr(eint1p(x, expx), z);
    1275        2779 :   set_avma(av); return z;
    1276             : }
    1277             : 
    1278             : static GEN
    1279         357 : cxeint1(GEN x, long prec)
    1280             : {
    1281         357 :   pari_sp av = avma, av2;
    1282             :   GEN q, S, run, z, H;
    1283         357 :   long n, E = prec2nbits(prec);
    1284             : 
    1285         357 :   if (gamma_use_asymp(x, E) && (z = eint1_asymp(x, NULL, prec))) return z;
    1286         252 :   E++;
    1287         252 :   if (gexpo(x) > 0)
    1288             :   { /* take cancellation into account, log2(\sum |x|^n / n!) = |x| / log(2) */
    1289          42 :     double dbx = dblmodulus(x);
    1290          42 :     long X = (long)((dbx + log(dbx))/M_LN2 + 10);
    1291          42 :     prec += nbits2extraprec(X);
    1292          42 :     x = gtofp(x, prec); E += X;
    1293             :   }
    1294         252 :   if (DEBUGLEVEL > 2) err_printf("eint1: using power series\n");
    1295         252 :   run = real_1(prec);
    1296         252 :   av2 = avma;
    1297         252 :   S = z = q = H = run;
    1298       48384 :   for (n = 2; gexpo(q) - gexpo(S) >= -E; n++)
    1299             :   {
    1300       48132 :     H = addrr(H, divru(run, n)); /* H = sum_{k<=n} 1/k */
    1301       48132 :     z = gdivgu(gmul(x,z), n);   /* z = x^(n-1)/n! */
    1302       48132 :     q = gmul(z, H); S = gadd(S, q);
    1303       48132 :     if ((n & 0x1ff) == 0) gerepileall(av2, 4, &z, &q, &S, &H);
    1304             :   }
    1305         252 :   S = gmul(gmul(x, S), gexp(gneg_i(x), prec));
    1306         252 :   return gerepileupto(av, gsub(S, gadd(glog(x, prec), mpeuler(prec))));
    1307             : }
    1308             : 
    1309             : GEN
    1310        3283 : eint1(GEN x, long prec)
    1311             : {
    1312        3283 :   switch(typ(x))
    1313             :   {
    1314         357 :     case t_COMPLEX: return cxeint1(x, prec);
    1315        2541 :     case t_REAL: break;
    1316         385 :     default: x = gtofp(x, prec);
    1317             :   }
    1318        2926 :   return mpeint1(x,NULL);
    1319             : }
    1320             : 
    1321             : GEN
    1322          49 : veceint1(GEN C, GEN nmax, long prec)
    1323             : {
    1324          49 :   if (!nmax) return eint1(C,prec);
    1325           7 :   if (typ(nmax) != t_INT) pari_err_TYPE("veceint1",nmax);
    1326           7 :   if (typ(C) != t_REAL) {
    1327           7 :     C = gtofp(C, prec);
    1328           7 :     if (typ(C) != t_REAL) pari_err_TYPE("veceint1",C);
    1329             :   }
    1330           7 :   if (signe(C) <= 0) pari_err_DOMAIN("veceint1", "argument", "<=", gen_0,C);
    1331           7 :   return mpveceint1(C, NULL, itos(nmax));
    1332             : }
    1333             : 
    1334             : /* j > 0, a t_REAL. Return sum_{m >= 0} a^m / j(j+1)...(j+m)).
    1335             :  * Stop when expo(summand) < E; note that s(j-1) = (a s(j) + 1) / (j-1). */
    1336             : static GEN
    1337         224 : mp_sum_j(GEN a, long j, long E, long prec)
    1338             : {
    1339         224 :   pari_sp av = avma;
    1340         224 :   GEN q = divru(real_1(prec), j), s = q;
    1341             :   long m;
    1342        4171 :   for (m = 0;; m++)
    1343             :   {
    1344        4171 :     if (expo(q) < E) break;
    1345        3947 :     q = mulrr(q, divru(a, m+j));
    1346        3947 :     s = addrr(s, q);
    1347             :   }
    1348         224 :   return gerepileuptoleaf(av, s);
    1349             : }
    1350             : /* Return the s_a(j), j <= J */
    1351             : static GEN
    1352         224 : sum_jall(GEN a, long J, long prec)
    1353             : {
    1354         224 :   GEN s = cgetg(J+1, t_VEC);
    1355         224 :   long j, E = -prec2nbits(prec) - 5;
    1356         224 :   gel(s, J) = mp_sum_j(a, J, E, prec);
    1357        9372 :   for (j = J-1; j; j--)
    1358        9148 :     gel(s,j) = divru(addrs(mulrr(a, gel(s,j+1)), 1), j);
    1359         224 :   return s;
    1360             : }
    1361             : 
    1362             : /* T a dense t_POL with t_REAL coeffs. Return T(n) [faster than poleval] */
    1363             : static GEN
    1364      361032 : rX_s_eval(GEN T, long n)
    1365             : {
    1366      361032 :   long i = lg(T)-1;
    1367      361032 :   GEN c = gel(T,i);
    1368     7466020 :   for (i--; i>=2; i--) c = gadd(mulrs(c,n),gel(T,i));
    1369      361032 :   return c;
    1370             : }
    1371             : 
    1372             : /* C>0 t_REAL, eC = exp(C). Return eint1(n*C) for 1<=n<=N. Absolute accuracy */
    1373             : GEN
    1374         231 : mpveceint1(GEN C, GEN eC, long N)
    1375             : {
    1376         231 :   const long prec = realprec(C);
    1377         231 :   long Nmin = 15; /* >= 1. E.g. between 10 and 30, but little effect */
    1378         231 :   GEN en, v, w = cgetg(N+1, t_VEC);
    1379             :   pari_sp av0;
    1380             :   double DL;
    1381             :   long n, j, jmax, jmin;
    1382         231 :   if (!N) return w;
    1383      364658 :   for (n = 1; n <= N; n++) gel(w,n) = cgetr(prec);
    1384         231 :   av0 = avma;
    1385         231 :   if (N < Nmin) Nmin = N;
    1386         231 :   if (!eC) eC = mpexp(C);
    1387         231 :   en = eC; affrr(eint1p(C, en), gel(w,1));
    1388        3395 :   for (n = 2; n <= Nmin; n++)
    1389             :   {
    1390             :     pari_sp av2;
    1391        3164 :     en = mulrr(en,eC); /* exp(n C) */
    1392        3164 :     av2 = avma;
    1393        3164 :     affrr(eint1p(mulru(C,n), en), gel(w,n));
    1394        3164 :     set_avma(av2);
    1395             :   }
    1396         231 :   if (Nmin == N) { set_avma(av0); return w; }
    1397             : 
    1398         224 :   DL = prec2nbits_mul(prec, M_LN2) + 5;
    1399         224 :   jmin = ceil(DL/log((double)N)) + 1;
    1400         224 :   jmax = ceil(DL/log((double)Nmin)) + 1;
    1401         224 :   v = sum_jall(C, jmax, prec);
    1402         224 :   en = powrs(eC, -N); /* exp(-N C) */
    1403         224 :   affrr(eint1p(mulru(C,N), invr(en)), gel(w,N));
    1404        5887 :   for (j = jmin, n = N-1; j <= jmax; j++)
    1405             :   {
    1406        5663 :     long limN = maxss((long)ceil(exp(DL/j)), Nmin);
    1407             :     GEN polsh;
    1408        5663 :     setlg(v, j+1);
    1409        5663 :     polsh = RgV_to_RgX_reverse(v, 0);
    1410      366695 :     for (; n >= limN; n--)
    1411             :     {
    1412      361032 :       pari_sp av2 = avma;
    1413      361032 :       GEN S = divri(mulrr(en, rX_s_eval(polsh, -n)), powuu(n,j));
    1414             :       /* w[n+1] - exp(-n C) * polsh(-n) / (-n)^j */
    1415      361032 :       GEN c = odd(j)? addrr(gel(w,n+1), S) : subrr(gel(w,n+1), S);
    1416      361032 :       affrr(c, gel(w,n)); set_avma(av2);
    1417      361032 :       en = mulrr(en,eC); /* exp(-n C) */
    1418             :     }
    1419             :   }
    1420         224 :   set_avma(av0); return w;
    1421             : }
    1422             : 
    1423             : /* erfc via numerical integration : assume real(x)>=1 */
    1424             : static GEN
    1425          14 : cxerfc_r1(GEN x, long prec)
    1426             : {
    1427             :   GEN h, h2, eh2, denom, res, lambda;
    1428             :   long u, v;
    1429          14 :   const double D = prec2nbits_mul(prec, M_LN2);
    1430          14 :   const long npoints = (long)ceil(D/M_PI)+1;
    1431          14 :   pari_sp av = avma;
    1432             :   {
    1433          14 :     double t = exp(-2*M_PI*M_PI/D); /* ~exp(-2*h^2) */
    1434          14 :     v = 30; /* bits that fit in both long and double mantissa */
    1435          14 :     u = (long)floor(t*(1L<<v));
    1436             :     /* define exp(-2*h^2) to be u*2^(-v) */
    1437             :   }
    1438          14 :   incrprec(prec);
    1439          14 :   x = gtofp(x,prec);
    1440          14 :   eh2 = sqrtr_abs(rtor(shiftr(dbltor(u),-v),prec));
    1441          14 :   h2 = negr(logr_abs(eh2));
    1442          14 :   h = sqrtr_abs(h2);
    1443          14 :   lambda = gdiv(x,h);
    1444          14 :   denom = gsqr(lambda);
    1445             :   { /* res = h/x + 2*x*h*sum(k=1,npoints,exp(-(k*h)^2)/(lambda^2+k^2)); */
    1446             :     GEN Uk; /* = exp(-(kh)^2) */
    1447          14 :     GEN Vk = eh2;/* = exp(-(2k+1)h^2) */
    1448          14 :     pari_sp av2 = avma;
    1449             :     long k;
    1450             :     /* k = 0 moved out for efficiency */
    1451          14 :     denom = gaddsg(1,denom);
    1452          14 :     Uk = Vk;
    1453          14 :     Vk = mulur(u,Vk); shiftr_inplace(Vk, -v);
    1454          14 :     res = gdiv(Uk, denom);
    1455         420 :     for (k = 1; k < npoints; k++)
    1456             :     {
    1457         406 :       if ((k & 255) == 0) gerepileall(av2,4,&denom,&Uk,&Vk,&res);
    1458         406 :       denom = gaddsg(2*k+1,denom);
    1459         406 :       Uk = mpmul(Uk,Vk);
    1460         406 :       Vk = mulur(u,Vk); shiftr_inplace(Vk, -v);
    1461         406 :       res = gadd(res, gdiv(Uk, denom));
    1462             :     }
    1463             :   }
    1464          14 :   res = gmul(res, gshift(lambda,1));
    1465             :   /* 0 term : */
    1466          14 :   res = gadd(res, ginv(lambda));
    1467          14 :   res = gmul(res, gdiv(gexp(gneg(gsqr(x)), prec), mppi(prec)));
    1468          14 :   if (rtodbl(real_i(x)) < sqrt(D))
    1469             :   {
    1470          14 :     GEN t = gmul(divrr(Pi2n(1,prec),h), x);
    1471          14 :     res = gsub(res, gdivsg(2, cxexpm1(t, prec)));
    1472             :   }
    1473          14 :   return gerepileupto(av,res);
    1474             : }
    1475             : 
    1476             : static GEN
    1477           7 : sererfc(GEN x, long prec)
    1478             : {
    1479           7 :   GEN u, z = invr(sqrtr_abs(Pi2n(-2,prec)));
    1480           7 :   setsigne(z, -1); /* -2/sqrt(Pi) */
    1481           7 :   z = gmul(z, integser(gmul(derivser(x), gexp(gneg(gsqr(x)), prec))));
    1482           7 :   u = polcoef_i(x, 0, varn(x));
    1483           7 :   if (!gcmp0(u)) z = gadd(z, gerfc(u, prec));
    1484           7 :   return z;
    1485             : }
    1486             : 
    1487             : GEN
    1488          70 : gerfc(GEN x, long prec)
    1489             : {
    1490             :   GEN z, xr, xi, res;
    1491             :   long s;
    1492             :   pari_sp av;
    1493             : 
    1494          70 :   switch(typ(x))
    1495             :   {
    1496          63 :     case t_INT: case t_REAL: case t_FRAC: case t_COMPLEX:
    1497          63 :       break;
    1498           7 :     default:
    1499           7 :       av = avma;
    1500           7 :       if ((z = toser_i(x))) return gerepileupto(av, sererfc(z,prec));
    1501           0 :       return trans_eval("erfc",gerfc,x,prec);
    1502             :   }
    1503             :   /* x a complex scalar */
    1504          63 :   x = trans_fix_arg(&prec,&x,&xr,&xi, &av,&res);
    1505          63 :   s = signe(xr);
    1506          63 :   if (s > 0 || (s == 0 && signe(xi) >= 0)) {
    1507          49 :     if (cmprs(xr, 1) > 0) /* use numerical integration */
    1508          14 :       z = cxerfc_r1(x, prec);
    1509             :     else
    1510             :     { /* erfc(x) = incgam(1/2,x^2)/sqrt(Pi) */
    1511          35 :       GEN sqrtpi = sqrtr(mppi(prec));
    1512          35 :       z = incgam0(ghalf, gsqr(x), sqrtpi, prec);
    1513          35 :       z = gdiv(z, sqrtpi);
    1514             :     }
    1515             :   }
    1516             :   else
    1517             :   { /* erfc(-x)=2-erfc(x) */
    1518             :     /* FIXME could decrease prec
    1519             :     long size = nbits2extraprec((imag(x)^2-real(x)^2)/log(2));
    1520             :     prec = size > 0 ? prec : prec + size;
    1521             :     */
    1522             :     /* NOT gsubsg(2, ...) : would create a result of
    1523             :      * huge accuracy if re(x)>>1, rounded to 2 by subsequent affc_fixlg... */
    1524          14 :     z = gsub(real2n(1,prec+EXTRAPREC64), gerfc(gneg(x), prec));
    1525             :   }
    1526          63 :   set_avma(av); return affc_fixlg(z, res);
    1527             : }
    1528             : 
    1529             : /***********************************************************************/
    1530             : /**                                                                   **/
    1531             : /**                      RIEMANN ZETA FUNCTION                        **/
    1532             : /**                                                                   **/
    1533             : /***********************************************************************/
    1534             : static const double log2PI = 1.83787706641;
    1535             : 
    1536             : static double
    1537        4578 : get_xinf(double beta)
    1538             : {
    1539        4578 :   const double maxbeta = 0.06415003; /* 3^(-2.5) */
    1540             :   double x0, y0, x1;
    1541             : 
    1542        4578 :   if (beta < maxbeta) return beta + pow(3*beta, 1.0/3.0);
    1543        4578 :   x0 = beta + M_PI/2.0;
    1544             :   for(;;)
    1545             :   {
    1546        7532 :     y0 = x0*x0;
    1547        7532 :     x1 = (beta+atan(x0)) * (1+y0) / y0 - 1/x0;
    1548        7532 :     if (0.99*x0 < x1) return x1;
    1549        2954 :     x0 = x1;
    1550             :   }
    1551             : }
    1552             : /* optimize for zeta( s + it, prec ), assume |s-1| > 0.1
    1553             :  * (if gexpo(u = s-1) < -5, we use the functional equation s->1-s) */
    1554             : static int
    1555        8610 : optim_zeta(GEN S, long prec, long *pp, long *pn)
    1556             : {
    1557             :   double s, t, alpha, beta, n, B;
    1558             :   long p;
    1559        8610 :   if (typ(S) == t_REAL) {
    1560         266 :     t = 0.;
    1561         266 :     s = rtodbl(S);
    1562             :   } else {
    1563        8344 :     t = fabs( rtodbl(gel(S,2)) );
    1564        8344 :     if (t > 2500) return 0; /* lfunlarge */
    1565        8309 :     s = rtodbl(gel(S,1));
    1566             :   }
    1567             : 
    1568        8575 :   B = prec2nbits_mul(prec, M_LN2);
    1569        8575 :   if (s > 0 && !t) /* positive real input */
    1570             :   {
    1571         238 :     beta = B + 0.61 + s*(log2PI - log(s));
    1572         238 :     if (beta > 0)
    1573             :     {
    1574         231 :       p = (long)ceil(beta / 2.0);
    1575         231 :       n = fabs(s + 2*p-1)/(2*M_PI);
    1576             :     }
    1577             :     else
    1578             :     {
    1579           7 :       p = 0;
    1580           7 :       n = exp((B - M_LN2) / s);
    1581             :     }
    1582             :   }
    1583        8337 :   else if (s <= 0 || t < 0.01) /* s < 0 may occur if s ~ 0 */
    1584        3752 :   { /* TODO: the crude bounds below are generally valid. Optimize ? */
    1585        3752 :     double l,l2, la = 1.; /* heuristic */
    1586        3752 :     double rlog, ilog; dcxlog(s-1,t, &rlog,&ilog);
    1587        3752 :     l2 = (s - 0.5)*rlog - t*ilog; /* = Re( (S - 1/2) log (S-1) ) */
    1588        3752 :     l = (B - l2 + s*log2PI) / (2. * (1.+ log((double)la)));
    1589        3752 :     l2 = dabs(s, t)/2;
    1590        3752 :     if (l < l2) l = l2;
    1591        3752 :     p = (long) ceil(l); if (p < 2) p = 2;
    1592        3752 :     n = 1 + dabs(p+s/2.-.25, t/2) * la / M_PI;
    1593             :   }
    1594             :   else
    1595             :   {
    1596        4585 :     double sn = dabs(s, t), L = log(sn/s);
    1597        4585 :     alpha = B - 0.39 + L + s*(log2PI - log(sn));
    1598        4585 :     beta = (alpha+s)/t - atan(s/t);
    1599        4585 :     p = 0;
    1600        4585 :     if (beta > 0)
    1601             :     {
    1602        4578 :       beta = 1.0 - s + t * get_xinf(beta);
    1603        4578 :       if (beta > 0) p = (long)ceil(beta / 2.0);
    1604             :     }
    1605             :     else
    1606           7 :       if (s < 1.0) p = 1;
    1607        4585 :     n = p? dabs(s + 2*p-1, t) / (2*M_PI) : exp((B-M_LN2+L) / s);
    1608             :   }
    1609        8575 :   *pp = p;
    1610        8575 :   *pn = (long)ceil(n);
    1611        8575 :   if (*pp < 0 || *pn < 0) pari_err_OVERFLOW("zeta");
    1612        8575 :   return 1;
    1613             : }
    1614             : 
    1615             : /* zeta(a*j+b), j=0..N-1, b>1, using sumalt. Johansonn's thesis, Algo 4.7.1 */
    1616             : static GEN
    1617         698 : veczetas(long a, long b, long N, long prec)
    1618             : {
    1619         698 :   const long n = ceil(2 + prec2nbits_mul(prec, M_LN2/1.7627));
    1620         698 :   pari_sp av = avma;
    1621         698 :   GEN c, d, z = zerovec(N);
    1622             :   long j, k;
    1623         698 :   c = d = int2n(2*n-1);
    1624       87066 :   for (k = n; k > 1; k--)
    1625             :   {
    1626       86368 :     GEN u, t = divii(d, powuu(k,b));
    1627       86361 :     if (!odd(k)) t = negi(t);
    1628       86393 :     gel(z,1) = addii(gel(z,1), t);
    1629       86426 :     u = powuu(k,a);
    1630     4097984 :     for (j = 1; j < N; j++)
    1631             :     {
    1632     4055749 :       t = divii(t,u); if (!signe(t)) break;
    1633     4008233 :       gel(z,j+1) = addii(gel(z,j+1), t);
    1634             :     }
    1635       87604 :     c = muluui(k,2*k-1,c);
    1636       86424 :     c = diviuuexact(c, 2*(n-k+1),n+k-1);
    1637       86385 :     d = addii(d,c);
    1638       86366 :     if (gc_needed(av,3))
    1639             :     {
    1640           7 :       if(DEBUGMEM>1) pari_warn(warnmem,"zetaBorwein, k = %ld", k);
    1641           7 :       gerepileall(av, 3, &c,&d,&z);
    1642             :     }
    1643             :   }
    1644             :   /* k = 1 */
    1645       53262 :   for (j = 1; j <= N; j++) gel(z,j) = addii(gel(z,j), d);
    1646         699 :   d = addiu(d, 1);
    1647       53320 :   for (j = 0, k = b - 1; j < N; j++, k += a)
    1648       52622 :     gel(z,j+1) = rdivii(shifti(gel(z,j+1), k), subii(shifti(d,k), d), prec);
    1649         698 :   return z;
    1650             : }
    1651             : /* zeta(a*j+b), j=0..N-1, b>1, using sumalt */
    1652             : GEN
    1653         217 : veczeta(GEN a, GEN b, long N, long prec)
    1654             : {
    1655         217 :   pari_sp av = avma;
    1656             :   long n, j, k;
    1657         217 :   GEN L, c, d, z = zerovec(N);
    1658         217 :   if (typ(a) == t_INT && typ(b) == t_INT)
    1659         126 :     return gerepilecopy(av, veczetas(itos(a),  itos(b), N, prec));
    1660          91 :   n = ceil(2 + prec2nbits_mul(prec, M_LN2/1.7627));
    1661          91 :   c = d = int2n(2*n-1);
    1662       13819 :   for (k = n; k; k--)
    1663             :   {
    1664             :     GEN u, t;
    1665       13728 :     L = logr_abs(utor(k, prec)); /* log(k) */
    1666       13728 :     t = gdiv(d, gexp(gmul(b, L), prec)); /* d / k^b */
    1667       13728 :     if (!odd(k)) t = gneg(t);
    1668       13728 :     gel(z,1) = gadd(gel(z,1), t);
    1669       13728 :     u = gexp(gmul(a, L), prec);
    1670      896888 :     for (j = 1; j < N; j++)
    1671             :     {
    1672      889583 :       t = gdiv(t,u); if (gexpo(t) < 0) break;
    1673      883160 :       gel(z,j+1) = gadd(gel(z,j+1), t);
    1674             :     }
    1675       13728 :     c = muluui(k,2*k-1,c);
    1676       13728 :     c = diviuuexact(c, 2*(n-k+1),n+k-1);
    1677       13728 :     d = addii(d,c);
    1678       13728 :     if (gc_needed(av,3))
    1679             :     {
    1680           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"veczeta, k = %ld", k);
    1681           0 :       gerepileall(av, 3, &c,&d,&z);
    1682             :     }
    1683             :   }
    1684          91 :   L = mplog2(prec);
    1685        5796 :   for (j = 0; j < N; j++)
    1686             :   {
    1687        5705 :     GEN u = gsubgs(gadd(b, gmulgu(a,j)), 1);
    1688        5705 :     GEN w = gexp(gmul(u, L), prec); /* 2^u */
    1689        5705 :     gel(z,j+1) = gdiv(gmul(gel(z,j+1), w), gmul(d,gsubgs(w,1)));
    1690             :   }
    1691          91 :   return gerepilecopy(av, z);
    1692             : }
    1693             : 
    1694             : GEN
    1695       27423 : constzeta(long n, long prec)
    1696             : {
    1697       27423 :   GEN o = zetazone, z;
    1698       27423 :   long l = o? lg(o): 0;
    1699             :   pari_sp av;
    1700       27423 :   if (l > n)
    1701             :   {
    1702       27022 :     long p = realprec(gel(o,1));
    1703       27022 :     if (p >= prec) return o;
    1704             :   }
    1705         572 :   n = maxss(n, l + 15);
    1706         572 :   av = avma; z = veczetas(1, 2, n-1, prec);
    1707         572 :   zetazone = gclone(vec_prepend(z, mpeuler(prec)));
    1708         572 :   set_avma(av); guncloneNULL(o); return zetazone;
    1709             : }
    1710             : 
    1711             : /* zeta(s) using sumalt, case h=0,N=1. Assume s > 1 */
    1712             : static GEN
    1713         598 : zetaBorwein(long s, long prec)
    1714             : {
    1715         598 :   pari_sp av = avma;
    1716         598 :   const long n = ceil(2 + prec2nbits_mul(prec, M_LN2/1.7627));
    1717             :   long k;
    1718         598 :   GEN c, d, z = gen_0;
    1719         598 :   c = d = int2n(2*n-1);
    1720      135672 :   for (k = n; k; k--)
    1721             :   {
    1722      135074 :     GEN t = divii(d, powuu(k,s));
    1723      135074 :     z = odd(k)? addii(z,t): subii(z,t);
    1724      135074 :     c = muluui(k,2*k-1,c);
    1725      135074 :     c = diviuuexact(c, 2*(n-k+1),n+k-1);
    1726      135074 :     d = addii(d,c);
    1727      135074 :     if (gc_needed(av,3))
    1728             :     {
    1729           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"zetaBorwein, k = %ld", k);
    1730           0 :       gerepileall(av, 3, &c,&d,&z);
    1731             :     }
    1732             :   }
    1733         598 :   return rdivii(shifti(z, s-1), subii(shifti(d,s-1), d), prec);
    1734             : }
    1735             : 
    1736             : /* assume k != 1 */
    1737             : GEN
    1738        4795 : szeta(long k, long prec)
    1739             : {
    1740        4795 :   pari_sp av = avma;
    1741             :   GEN z;
    1742             : 
    1743        4795 :   if (!k) { z = real2n(-1, prec); setsigne(z,-1); return z; }
    1744        4788 :   if (k < 0)
    1745             :   {
    1746         252 :     if (!odd(k)) return gen_0;
    1747             :     /* the one value such that k < 0 and 1 - k < 0, due to overflow */
    1748         252 :     if ((ulong)k == (HIGHBIT | 1))
    1749           0 :       pari_err_OVERFLOW("zeta [large negative argument]");
    1750         252 :     k = 1-k;
    1751         252 :     z = bernreal(k, prec); togglesign(z);
    1752         252 :     return gerepileuptoleaf(av, divru(z, k));
    1753             :   }
    1754             :   /* k > 1 */
    1755        4536 :   if (k > prec2nbits(prec)+1) return real_1(prec);
    1756        4529 :   if (zetazone && realprec(gel(zetazone,1)) >= prec && lg(zetazone) > k)
    1757         119 :     return rtor(gel(zetazone, k), prec);
    1758        4410 :   if (!odd(k))
    1759             :   {
    1760             :     GEN B;
    1761        2639 :     if (!bernzone) constbern(0);
    1762        2639 :     if (k < lg(bernzone))
    1763        2184 :       B = gel(bernzone, k>>1);
    1764             :     else
    1765             :     {
    1766         455 :       if (bernbitprec(k) > prec2nbits(prec))
    1767           0 :         return gerepileupto(av, invr(inv_szeta_euler(k, prec)));
    1768         455 :       B = bernfrac(k);
    1769             :     }
    1770             :     /* B = B_k */
    1771        2639 :     z = gmul(powru(Pi2n(1, prec + EXTRAPREC64), k), B);
    1772        2639 :     z = divrr(z, mpfactr(k,prec));
    1773        2639 :     setsigne(z, 1); shiftr_inplace(z, -1);
    1774             :   }
    1775             :   else
    1776             :   {
    1777        1771 :     double p = prec2nbits_mul(prec,0.393); /* bit / log_2(3+sqrt(8)) */
    1778        1771 :     p = log2(p * log(p));
    1779        3542 :     z = (p * k > prec2nbits(prec))? invr(inv_szeta_euler(k, prec))
    1780        1771 :                                   : zetaBorwein(k, prec);
    1781             :   }
    1782        4410 :   return gerepileuptoleaf(av, z);
    1783             : }
    1784             : 
    1785             : /* Ensure |1-s| >= 1/32 and (|s| <= 1/32 or real(s) >= 1/2) */
    1786             : static int
    1787        8624 : zeta_funeq(GEN *ps)
    1788             : {
    1789        8624 :   GEN s = *ps, u;
    1790        8624 :   if (typ(s) == t_REAL)
    1791             :   {
    1792         273 :     u = subsr(1, s);
    1793         273 :     if (expo(u) >= -5
    1794         266 :         && ((signe(s) > 0 && expo(s) >= -1) || expo(s) <= -5)) return 0;
    1795             :   }
    1796             :   else
    1797             :   {
    1798        8351 :     GEN sig = gel(s,1);
    1799        8351 :     if (fabs(rtodbl(gel(s,2))) > 2500) return 0; /* lfunlarge */
    1800        8316 :     u = gsubsg(1, s);
    1801        8316 :     if (gexpo(u) >= -5
    1802        8309 :         && ((signe(sig) > 0 && expo(sig) >= -1) || gexpo(s) <= -5)) return 0;
    1803             :   }
    1804        1505 :   *ps = u; return 1;
    1805             : }
    1806             : /* s0 a t_INT, t_REAL or t_COMPLEX.
    1807             :  * If a t_INT, assume it's not a trivial case (i.e we have s0 > 1, odd) */
    1808             : static GEN
    1809        8631 : czeta(GEN s0, long prec)
    1810             : {
    1811        8631 :   GEN ms, s, u, y, res, tes, sig, tau, invn2, ns, Ns, funeq_factor = NULL;
    1812             :   long i, nn, lim, lim2;
    1813        8631 :   pari_sp av0 = avma, av, av2;
    1814             :   pari_timer T;
    1815             : 
    1816        8631 :   if (DEBUGLEVEL>2) timer_start(&T);
    1817        8631 :   s = trans_fix_arg(&prec,&s0,&sig,&tau,&av,&res);
    1818        8631 :   if (typ(s0) == t_INT) return gerepileupto(av0, gzeta(s0, prec));
    1819        8624 :   if (zeta_funeq(&s)) /* s -> 1-s */
    1820             :   { /* Gamma(s) (2Pi)^-s 2 cos(Pi s/2) [new s] */
    1821        1505 :     GEN t = gmul(ggamma(s,prec), pow2Pis(gsubgs(s0,1), prec));
    1822        1505 :     sig = real_i(s);
    1823        1505 :     funeq_factor = gmul2n(gmul(t, gsin(gmul(Pi2n(-1,prec),s0), prec)), 1);
    1824             :   }
    1825        8624 :   if (gcmpgs(sig, prec2nbits(prec) + 1) > 0) { /* zeta(s) = 1 */
    1826          14 :     if (!funeq_factor) { set_avma(av0); return real_1(prec); }
    1827           7 :     return gerepileupto(av0, funeq_factor);
    1828             :   }
    1829        8610 :   if (!optim_zeta(s, prec, &lim, &nn))
    1830             :   {
    1831          35 :     long bit = prec2nbits(prec);
    1832          35 :     y = lfun(lfuninit(gen_1, cgetg(1,t_VEC), 0, bit), s, bit);
    1833          35 :     if (funeq_factor) y = gmul(y, funeq_factor);
    1834          35 :     set_avma(av); return affc_fixlg(y,res);
    1835             :   }
    1836        8575 :   if (DEBUGLEVEL>2) err_printf("lim, nn: [%ld, %ld]\n", lim, nn);
    1837        8575 :   ms = gneg(s);
    1838        8575 :   if (umuluu_le(nn, prec, 10000000))
    1839             :   {
    1840        8575 :     incrprec(prec); /* one extra word of precision */
    1841        8575 :     Ns = vecpowug(nn, ms, prec);
    1842        8575 :     ns = gel(Ns,nn); setlg(Ns, nn);
    1843        8575 :     y = gadd(gmul2n(ns, -1), RgV_sum(Ns));
    1844             :   }
    1845             :   else
    1846             :   {
    1847           0 :     Ns = dirpowerssum(nn, ms, 0, prec);
    1848           0 :     incrprec(prec); /* one extra word of precision */
    1849           0 :     ns = gpow(utor(nn, prec), ms, prec);
    1850           0 :     y = gsub(Ns, gmul2n(ns, -1));
    1851             :   }
    1852        8575 :   if (DEBUGLEVEL>2) timer_printf(&T,"sum from 1 to N");
    1853        8575 :   constbern(lim);
    1854        8575 :   if (DEBUGLEVEL>2) timer_start(&T);
    1855        8575 :   invn2 = divri(real_1(prec), sqru(nn)); lim2 = lim<<1;
    1856        8575 :   tes = bernfrac(lim2);
    1857             :   {
    1858             :     GEN s1, s2, s3, s4, s5;
    1859        8575 :     s2 = gmul(s, gsubgs(s,1));
    1860        8575 :     s3 = gmul2n(invn2,3);
    1861        8575 :     av2 = avma;
    1862        8575 :     s1 = gsubgs(gmul2n(s,1), 1);
    1863        8575 :     s4 = gmul(invn2, gmul2n(gaddsg(4*lim-2,s1),1));
    1864        8575 :     s5 = gmul(invn2, gadd(s2, gmulsg(lim2, gaddgs(s1, lim2))));
    1865      805263 :     for (i = lim2-2; i>=2; i -= 2)
    1866             :     {
    1867      796688 :       s5 = gsub(s5, s4);
    1868      796688 :       s4 = gsub(s4, s3);
    1869      796688 :       tes = gadd(bernfrac(i), gdivgunextu(gmul(s5,tes), i+1));
    1870      796688 :       if (gc_needed(av2,3))
    1871             :       {
    1872           0 :         if(DEBUGMEM>1) pari_warn(warnmem,"czeta i = %ld", i);
    1873           0 :         gerepileall(av2,3, &tes,&s5,&s4);
    1874             :       }
    1875             :     }
    1876        8575 :     u = gmul(gmul(tes,invn2), gmul2n(s2, -1));
    1877        8575 :     tes = gmulsg(nn, gaddsg(1, u));
    1878             :   }
    1879        8575 :   if (DEBUGLEVEL>2) timer_printf(&T,"Bernoulli sum");
    1880             :   /* y += tes n^(-s) / (s-1) */
    1881        8575 :   y = gadd(y, gmul(tes, gdiv(ns, gsubgs(s,1))));
    1882        8575 :   if (funeq_factor) y = gmul(y, funeq_factor);
    1883        8575 :   set_avma(av); return affc_fixlg(y,res);
    1884             : }
    1885             : /* v a t_VEC/t_COL */
    1886             : int
    1887          14 : RgV_is_arithprog(GEN v, GEN *a, GEN *b)
    1888             : {
    1889          14 :   pari_sp av = avma, av2;
    1890          14 :   long i, n = lg(v)-1;
    1891          14 :   if (n == 0) { *a = *b = gen_0; return 1; }
    1892          14 :   *a = gel(v,1);
    1893          14 :   if (n == 1) { * b = gen_0; return 1; }
    1894          14 :   *b = gsub(gel(v,2), *a); av2 = avma;
    1895          28 :   for (i = 2; i < n; i++)
    1896          14 :     if (!gequal(*b, gsub(gel(v,i+1), gel(v,i)))) return gc_int(av,0);
    1897          14 :   return gc_int(av2,1);
    1898             : }
    1899             : 
    1900             : GEN
    1901       13440 : gzeta(GEN x, long prec)
    1902             : {
    1903       13440 :   pari_sp av = avma;
    1904             :   GEN y;
    1905       13440 :   if (gequal1(x)) pari_err_DOMAIN("zeta", "argument", "=", gen_1, x);
    1906       13419 :   switch(typ(x))
    1907             :   {
    1908        4557 :     case t_INT:
    1909        4557 :       if (is_bigint(x))
    1910             :       {
    1911          21 :         if (signe(x) > 0) return real_1(prec);
    1912          14 :         if (mod2(x) == 0) return real_0(prec);
    1913           7 :         pari_err_OVERFLOW("zeta [large negative argument]");
    1914             :       }
    1915        4536 :       return szeta(itos(x),prec);
    1916        8631 :     case t_REAL: case t_COMPLEX: return czeta(x,prec);
    1917          14 :     case t_PADIC: return Qp_zeta(x);
    1918          21 :     case t_VEC: case t_COL:
    1919             :     {
    1920             :       GEN a, b;
    1921          21 :       long n = lg(x) - 1;
    1922          21 :       if (n > 1 && RgV_is_arithprog(x, &a, &b))
    1923             :       {
    1924          14 :         if (!is_real_t(typ(a)) || !is_real_t(typ(b)) || gcmpgs(a, 1) <= 0)
    1925           0 :         { set_avma(av); break; }
    1926          14 :         a = veczeta(b, a, n, prec);
    1927          14 :         settyp(a, typ(x)); return a;
    1928             :       }
    1929             :     }
    1930             :     default:
    1931         203 :       if (!(y = toser_i(x))) break;
    1932          35 :       if (gequal1(y))
    1933           7 :         pari_err_DOMAIN("zeta", "argument", "=", gen_1, y);
    1934          28 :       return gerepileupto(av, lfun(gen_1,y,prec2nbits(prec)));
    1935             :   }
    1936         168 :   return trans_eval("zeta",gzeta,x,prec);
    1937             : }
    1938             : 
    1939             : /***********************************************************************/
    1940             : /**                                                                   **/
    1941             : /**                    FONCTIONS POLYLOGARITHME                       **/
    1942             : /**                                                                   **/
    1943             : /***********************************************************************/
    1944             : 
    1945             : /* smallish k such that bernbitprec(K) > bit + Kdz, K = 2k+4 */
    1946             : static long
    1947          21 : get_k(double dz, long bit)
    1948             : {
    1949             :   long a, b;
    1950          21 :   for (b = 128;; b <<= 1)
    1951          21 :     if (bernbitprec(b) > bit + b*dz) break;
    1952          21 :   if (b == 128) return 128;
    1953           0 :   a = b >> 1;
    1954           0 :   while (b - a > 64)
    1955             :   {
    1956           0 :     long c = (a+b) >> 1;
    1957           0 :     if (bernbitprec(c) > bit + c*dz) b = c; else a = c;
    1958             :   }
    1959           0 :   return b >> 1;
    1960             : }
    1961             : 
    1962             : /* m >= 2. Validity domain |log x| < 2*Pi, contains log |x| < 5.44,
    1963             :  * Li_m(x = e^z) = sum_{n >= 0} zeta(m-n) z^n / n!
    1964             :  *    with zeta(1) := H_{m-1} - log(-z) */
    1965             : static GEN
    1966          21 : cxpolylog(long m, GEN x, long prec)
    1967             : {
    1968             :   long li, n, k, ksmall, real;
    1969             :   GEN vz, z, Z, h, q, s, S;
    1970             :   pari_sp av;
    1971             :   double dz;
    1972             :   pari_timer T;
    1973             : 
    1974          21 :   if (gequal1(x)) return szeta(m,prec);
    1975             :   /* x real <= 1 ==> Li_m(x) real */
    1976          21 :   real = (typ(x) == t_REAL && (expo(x) < 0 || signe(x) <= 0));
    1977             : 
    1978          21 :   vz = constzeta(m, prec);
    1979          21 :   z = glog(x,prec);
    1980             :   /* n = 0 */
    1981          21 :   q = gen_1; s = gel(vz, m);
    1982          28 :   for (n=1; n < m-1; n++)
    1983             :   {
    1984           7 :     q = gdivgu(gmul(q,z),n);
    1985           7 :     s = gadd(s, gmul(gel(vz,m-n), real? real_i(q): q));
    1986             :   }
    1987             :   /* n = m-1 */
    1988          21 :     q = gdivgu(gmul(q,z),n); /* multiply by "zeta(1)" */
    1989          21 :     h = gmul(q, gsub(harmonic(m-1), glog(gneg_i(z),prec)));
    1990          21 :     s = gadd(s, real? real_i(h): h);
    1991             :   /* n = m */
    1992          21 :     q = gdivgu(gmul(q,z),m);
    1993          21 :     s = gadd(s, gdivgs(real? real_i(q): q, -2)); /* zeta(0) = -1/2 */
    1994             :   /* n = m+1 */
    1995          21 :     q = gdivgu(gmul(q,z),m+1); /* = z^(m+1) / (m+1)! */
    1996          21 :     s = gadd(s, gdivgs(real? real_i(q): q, -12)); /* zeta(-1) = -1/12 */
    1997             : 
    1998          21 :   li = -(prec2nbits(prec)+1);
    1999          21 :   if (DEBUGLEVEL) timer_start(&T);
    2000          21 :   dz = dbllog2(z) - log2PI; /*  ~ log2(|z|/2Pi) */
    2001             :   /* sum_{k >= 1} zeta(-1-2k) * z^(2k+m+1) / (2k+m+1)!
    2002             :    * = 2 z^(m-1) sum_{k >= 1} zeta(2k+2) * Z^(k+1) / (2k+2)..(2k+1+m)), where
    2003             :    * Z = -(z/2Pi)^2. Stop at 2k = (li - (m-1)*Lz - m) / dz, Lz = log2 |z| */
    2004             :   /* We cut the sum in two: small values of k first */
    2005          21 :   Z = gsqr(z); av = avma;
    2006          21 :   ksmall = get_k(dz, prec2nbits(prec));
    2007          21 :   constbern(ksmall);
    2008         469 :   for(k = 1; k < ksmall; k++)
    2009             :   {
    2010         469 :     GEN t = q = gdivgunextu(gmul(q,Z), 2*k+m); /* z^(2k+m+1)/(2k+m+1)! */
    2011         469 :     if (real) t = real_i(t);
    2012         469 :     t = gmul(t, gdivgu(bernfrac(2*k+2), 2*k+2)); /* - t * zeta(1-(2k+2)) */
    2013         469 :     s = gsub(s, t);
    2014         469 :     if (gexpo(t)  < li) return s;
    2015             :     /* large values ? */
    2016         448 :     if ((k & 0x1ff) == 0) gerepileall(av, 2, &s, &q);
    2017             :   }
    2018           0 :   if (DEBUGLEVEL>2) timer_printf(&T, "polylog: small k <= %ld", k);
    2019           0 :   Z = gneg(gsqr(gdiv(z, Pi2n(1,prec))));
    2020           0 :   q = gmul(gpowgs(z, m-1), gpowgs(Z, k+1)); /* Z^(k+1) * z^(m-1) */
    2021           0 :   S = gen_0; av = avma;
    2022           0 :   for(;; k++)
    2023           0 :   {
    2024           0 :     GEN t = q;
    2025             :     long b;
    2026           0 :     if (real) t = real_i(t);
    2027           0 :     b = prec + gexpo(t) / BITS_IN_LONG; /* decrease accuracy */
    2028           0 :     if (b == 2) break;
    2029             :     /* t * zeta(2k+2) / (2k+2)..(2k+1+m) */
    2030           0 :     t = gdiv(t, mulri(inv_szeta_euler(2*k+2, b),
    2031           0 :                       mulu_interval(2*k+2, 2*k+1+m)));
    2032           0 :     S = gadd(S, t); if (gexpo(t)  < li) break;
    2033           0 :     q = gmul(q, Z);
    2034           0 :     if ((k & 0x1ff) == 0) gerepileall(av, 2, &S, &q);
    2035             :   }
    2036           0 :   if (DEBUGLEVEL>2) timer_printf(&T, "polylog: large k <= %ld", k);
    2037           0 :   return gadd(s, gmul2n(S,1));
    2038             : }
    2039             : 
    2040             : static GEN
    2041          42 : Li1(GEN x, long prec) { return gneg(glog(gsubsg(1, x), prec)); }
    2042             : 
    2043             : static GEN
    2044         203 : polylog(long m, GEN x, long prec)
    2045             : {
    2046             :   long l, e, i, G, sx;
    2047             :   pari_sp av, av1;
    2048             :   GEN X, Xn, z, p1, p2, y, res;
    2049             : 
    2050         203 :   if (m < 0) pari_err_DOMAIN("polylog", "index", "<", gen_0, stoi(m));
    2051         203 :   if (!m) return mkfrac(gen_m1,gen_2);
    2052         203 :   if (gequal0(x)) return gcopy(x);
    2053         203 :   if (m==1) { av = avma; return gerepileupto(av, Li1(x, prec)); }
    2054             : 
    2055         168 :   l = precision(x);
    2056         168 :   if (!l) l = prec; else prec = l;
    2057         168 :   res = cgetc(l); av = avma;
    2058         168 :   x = gtofp(x, l+EXTRAPREC64);
    2059         168 :   e = gexpo(gnorm(x));
    2060         168 :   if (!e || e == -1) {
    2061          21 :     y = cxpolylog(m,x,prec);
    2062          21 :     set_avma(av); return affc_fixlg(y, res);
    2063             :   }
    2064         147 :   X = (e > 0)? ginv(x): x;
    2065         147 :   G = -prec2nbits(l);
    2066         147 :   av1 = avma;
    2067         147 :   y = Xn = X;
    2068       68159 :   for (i=2; ; i++)
    2069             :   {
    2070       68159 :     Xn = gmul(X,Xn); p2 = gdiv(Xn,powuu(i,m));
    2071       68159 :     y = gadd(y,p2);
    2072       68159 :     if (gexpo(p2) <= G) break;
    2073             : 
    2074       68012 :     if (gc_needed(av1,1))
    2075             :     {
    2076           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"polylog");
    2077           0 :       gerepileall(av1,2, &y, &Xn);
    2078             :     }
    2079             :   }
    2080         147 :   if (e < 0) { set_avma(av); return affc_fixlg(y, res); }
    2081             : 
    2082          28 :   sx = gsigne(imag_i(x));
    2083          28 :   if (!sx)
    2084             :   {
    2085          28 :     if (m&1) sx = gsigne(gsub(gen_1, real_i(x)));
    2086          21 :     else     sx = - gsigne(real_i(x));
    2087             :   }
    2088          28 :   z = divri(mppi(l), mpfact(m-1)); setsigne(z, sx);
    2089          28 :   z = mkcomplex(gen_0, z);
    2090             : 
    2091          28 :   if (m == 2)
    2092             :   { /* same formula as below, written more efficiently */
    2093          21 :     y = gneg_i(y);
    2094          21 :     if (typ(x) == t_REAL && signe(x) < 0)
    2095           7 :       p1 = logr_abs(x);
    2096             :     else
    2097          14 :       p1 = gsub(glog(x,l), z);
    2098          21 :     p1 = gmul2n(gsqr(p1), -1); /* = (log(-x))^2 / 2 */
    2099             : 
    2100          21 :     p1 = gadd(p1, divru(sqrr(mppi(l)), 6));
    2101          21 :     p1 = gneg_i(p1);
    2102             :   }
    2103             :   else
    2104             :   {
    2105           7 :     GEN logx = glog(x,l), logx2 = gsqr(logx), vz = constzeta(m, l);
    2106           7 :     p1 = mkfrac(gen_m1,gen_2);
    2107          14 :     for (i = m-2; i >= 0; i -= 2)
    2108           7 :       p1 = gadd(gel(vz, m-i), gmul(p1, gdivgunextu(logx2, i+1)));
    2109           7 :     if (m&1) p1 = gmul(logx,p1); else y = gneg_i(y);
    2110           7 :     p1 = gadd(gmul2n(p1,1), gmul(z,gpowgs(logx,m-1)));
    2111           7 :     if (typ(x) == t_REAL && signe(x) < 0) p1 = real_i(p1);
    2112             :   }
    2113          28 :   y = gadd(y,p1);
    2114          28 :   set_avma(av); return affc_fixlg(y, res);
    2115             : }
    2116             : static GEN
    2117         119 : RIpolylog(long m, GEN x, long real, long prec)
    2118             : {
    2119         119 :   GEN y = polylog(m, x, prec);
    2120         119 :   return real? real_i(y): imag_i(y);
    2121             : }
    2122             : GEN
    2123          21 : dilog(GEN x, long prec) { return gpolylog(2, x, prec); }
    2124             : 
    2125             : /* x a floating point number, t_REAL or t_COMPLEX of t_REAL */
    2126             : static GEN
    2127          42 : logabs(GEN x)
    2128             : {
    2129             :   GEN y;
    2130          42 :   if (typ(x) == t_COMPLEX)
    2131             :   {
    2132           7 :     y = logr_abs( cxnorm(x) );
    2133           7 :     shiftr_inplace(y, -1);
    2134             :   } else
    2135          35 :     y = logr_abs(x);
    2136          42 :   return y;
    2137             : }
    2138             : 
    2139             : static GEN
    2140          21 : polylogD(long m, GEN x, long flag, long prec)
    2141             : {
    2142          21 :   long fl = 0, k, l, m2;
    2143             :   pari_sp av;
    2144             :   GEN p1, p2, y;
    2145             : 
    2146          21 :   if (gequal0(x)) return gcopy(x);
    2147          21 :   m2 = m&1;
    2148          21 :   if (gequal1(x) && m>=2) return m2? szeta(m,prec): gen_0;
    2149          21 :   av = avma; l = precision(x);
    2150          21 :   if (!l) { l = prec; x = gtofp(x,l); }
    2151          21 :   p1 = logabs(x);
    2152          21 :   if (signe(p1) > 0) { x = ginv(x); fl = !m2; } else setabssign(p1);
    2153             :   /* |x| <= 1, p1 = - log|x| >= 0 */
    2154          21 :   p2 = gen_1;
    2155          21 :   y = RIpolylog(m, x, m2, l);
    2156          84 :   for (k = 1; k < m; k++)
    2157             :   {
    2158          63 :     GEN t = RIpolylog(m-k, x, m2, l);
    2159          63 :     p2 = gdivgu(gmul(p2,p1), k); /* (-log|x|)^k / k! */
    2160          63 :     y = gadd(y, gmul(p2, t));
    2161             :   }
    2162          21 :   if (m2)
    2163             :   {
    2164          14 :     p1 = flag? gdivgs(p1, -2*m): gdivgs(logabs(gsubsg(1,x)), m);
    2165          14 :     y = gadd(y, gmul(p2, p1));
    2166             :   }
    2167          21 :   if (fl) y = gneg(y);
    2168          21 :   return gerepileupto(av, y);
    2169             : }
    2170             : 
    2171             : static GEN
    2172          14 : polylogP(long m, GEN x, long prec)
    2173             : {
    2174          14 :   long fl = 0, k, l, m2;
    2175             :   pari_sp av;
    2176             :   GEN p1,y;
    2177             : 
    2178          14 :   if (gequal0(x)) return gcopy(x);
    2179          14 :   m2 = m&1;
    2180          14 :   if (gequal1(x) && m>=2) return m2? szeta(m,prec): gen_0;
    2181          14 :   av = avma; l = precision(x);
    2182          14 :   if (!l) { l = prec; x = gtofp(x,l); }
    2183          14 :   p1 = logabs(x);
    2184          14 :   if (signe(p1) > 0) { x = ginv(x); fl = !m2; setsigne(p1, -1); }
    2185             :   /* |x| <= 1 */
    2186          14 :   y = RIpolylog(m, x, m2, l);
    2187          14 :   if (m==1)
    2188             :   {
    2189           7 :     shiftr_inplace(p1, -1); /* log |x| / 2 */
    2190           7 :     y = gadd(y, p1);
    2191             :   }
    2192             :   else
    2193             :   { /* m >= 2, \sum_{0 <= k <= m} 2^k B_k/k! (log |x|)^k Li_{m-k}(x),
    2194             :        with Li_0(x) := -1/2 */
    2195           7 :     GEN u, t = RIpolylog(m-1, x, m2, l);
    2196           7 :     u = gneg_i(p1); /* u = 2 B1 log |x| */
    2197           7 :     y = gadd(y, gmul(u, t));
    2198           7 :     if (m > 2)
    2199             :     {
    2200             :       GEN p2;
    2201           7 :       shiftr_inplace(p1, 1); /* 2log|x| <= 0 */
    2202           7 :       constbern(m>>1);
    2203           7 :       p1 = sqrr(p1);
    2204           7 :       p2 = shiftr(p1,-1);
    2205          21 :       for (k = 2; k < m; k += 2)
    2206             :       {
    2207          14 :         if (k > 2) p2 = gdivgunextu(gmul(p2,p1),k-1); /* 2^k/k! log^k |x|*/
    2208          14 :         t = RIpolylog(m-k, x, m2, l);
    2209          14 :         u = gmul(p2, bernfrac(k));
    2210          14 :         y = gadd(y, gmul(u, t));
    2211             :       }
    2212             :     }
    2213             :   }
    2214          14 :   if (fl) y = gneg(y);
    2215          14 :   return gerepileupto(av, y);
    2216             : }
    2217             : 
    2218             : static GEN
    2219         175 : gpolylog_i(void *E, GEN x, long prec)
    2220             : {
    2221         175 :   pari_sp av = avma;
    2222         175 :   long i, n, v, m = (long)E;
    2223             :   GEN a, y;
    2224             : 
    2225         175 :   if (m <= 0)
    2226             :   {
    2227          28 :     a = gmul(x, poleval(eulerianpol(-m, 0), x));
    2228          28 :     return gerepileupto(av, gdiv(a, gpowgs(gsubsg(1, x), 1-m)));
    2229             :   }
    2230         147 :   switch(typ(x))
    2231             :   {
    2232          84 :     case t_REAL: case t_COMPLEX: return polylog(m,x,prec);
    2233           7 :     case t_INTMOD: case t_PADIC: pari_err_IMPL( "padic polylogarithm");
    2234          56 :     default:
    2235          56 :       av = avma; if (!(y = toser_i(x))) break;
    2236          21 :       if (!m) { set_avma(av); return mkfrac(gen_m1,gen_2); }
    2237          21 :       if (m==1) return gerepileupto(av, Li1(y, prec));
    2238          21 :       if (gequal0(y)) return gerepilecopy(av, y);
    2239          21 :       v = valser(y);
    2240          21 :       if (v < 0) pari_err_DOMAIN("polylog","valuation", "<", gen_0, x);
    2241          14 :       if (v > 0) {
    2242           7 :         n = (lg(y)-3 + v) / v;
    2243           7 :         a = zeroser(varn(y), lg(y)-2);
    2244          35 :         for (i=n; i>=1; i--)
    2245          28 :           a = gmul(y, gadd(a, powis(utoipos(i),-m)));
    2246             :       } else { /* v == 0 */
    2247           7 :         long vy = varn(y);
    2248           7 :         GEN a0 = polcoef_i(y, 0, -1), t = gdiv(derivser(y), y);
    2249           7 :         a = Li1(y, prec);
    2250          14 :         for (i=2; i<=m; i++)
    2251           7 :           a = gadd(gpolylog(i, a0, prec), integ(gmul(t, a), vy));
    2252             :       }
    2253          14 :       return gerepileupto(av, a);
    2254             :   }
    2255          35 :   return trans_evalgen("polylog", E, gpolylog_i, x, prec);
    2256             : }
    2257             : GEN
    2258         133 : gpolylog(long m, GEN x, long prec) { return gpolylog_i((void*)m, x, prec); }
    2259             : 
    2260             : GEN
    2261         147 : polylog0(long m, GEN x, long flag, long prec)
    2262             : {
    2263         147 :   switch(flag)
    2264             :   {
    2265         105 :     case 0: return gpolylog(m,x,prec);
    2266          14 :     case 1: return polylogD(m,x,0,prec);
    2267           7 :     case 2: return polylogD(m,x,1,prec);
    2268          14 :     case 3: return polylogP(m,x,prec);
    2269           7 :     default: pari_err_FLAG("polylog");
    2270             :   }
    2271             :   return NULL; /* LCOV_EXCL_LINE */
    2272             : }
    2273             : 
    2274             : GEN
    2275      108511 : upper_to_cx(GEN x, long *prec)
    2276             : {
    2277      108511 :   long tx = typ(x), l;
    2278      108511 :   if (tx == t_QUAD) { x = quadtofp(x, *prec); tx = typ(x); }
    2279      108511 :   switch(tx)
    2280             :   {
    2281      108490 :     case t_COMPLEX:
    2282      108490 :       if (gsigne(gel(x,2)) > 0) break; /*fall through*/
    2283             :     case t_REAL: case t_INT: case t_FRAC:
    2284          14 :       pari_err_DOMAIN("modular function", "Im(argument)", "<=", gen_0, x);
    2285           7 :     default:
    2286           7 :       pari_err_TYPE("modular function", x);
    2287             :   }
    2288      108490 :   l = precision(x); if (l) *prec = l;
    2289      108490 :   return x;
    2290             : }
    2291             : 
    2292             : /* sqrt(3)/2 */
    2293             : static GEN
    2294        2233 : sqrt32(long prec) { GEN z = sqrtr_abs(utor(3,prec)); setexpo(z, -1); return z; }
    2295             : /* exp(i x), x = k pi/12 */
    2296             : static GEN
    2297        3657 : e12(ulong k, long prec)
    2298             : {
    2299             :   int s, sPi, sPiov2;
    2300             :   GEN z, t;
    2301        3657 :   k %= 24;
    2302        3657 :   if (!k) return gen_1;
    2303        3657 :   if (k == 12) return gen_m1;
    2304        3657 :   if (k >12) { s = 1; k = 24 - k; } else s = 0; /* x -> 2pi - x */
    2305        3657 :   if (k > 6) { sPi = 1; k = 12 - k; } else sPi = 0; /* x -> pi  - x */
    2306        3657 :   if (k > 3) { sPiov2 = 1; k = 6 - k; } else sPiov2 = 0; /* x -> pi/2 - x */
    2307        3657 :   z = cgetg(3, t_COMPLEX);
    2308        3657 :   switch(k)
    2309             :   {
    2310         716 :     case 0: gel(z,1) = icopy(gen_1); gel(z,2) = gen_0; break;
    2311         812 :     case 1: t = gmul2n(addrs(sqrt32(prec), 1), -1);
    2312         812 :       gel(z,1) = sqrtr(t);
    2313         812 :       gel(z,2) = gmul2n(invr(gel(z,1)), -2); break;
    2314             : 
    2315        1421 :     case 2: gel(z,1) = sqrt32(prec);
    2316        1421 :             gel(z,2) = real2n(-1, prec); break;
    2317             : 
    2318         708 :     case 3: gel(z,1) = sqrtr_abs(real2n(-1,prec));
    2319         708 :             gel(z,2) = rcopy(gel(z,1)); break;
    2320             :   }
    2321        3657 :   if (sPiov2) swap(gel(z,1), gel(z,2));
    2322        3657 :   if (sPi) togglesign(gel(z,1));
    2323        3657 :   if (s)   togglesign(gel(z,2));
    2324        3657 :   return z;
    2325             : }
    2326             : /* z a t_FRAC */
    2327             : static GEN
    2328       15052 : expIPifrac(GEN z, long prec)
    2329             : {
    2330       15052 :   GEN n = gel(z,1), d = gel(z,2);
    2331       15052 :   ulong r, q = uabsdivui_rem(12, d, &r);
    2332       15052 :   if (!r) return e12(q * umodiu(n, 24), prec); /* 12 | d */
    2333       11395 :   n = centermodii(n, shifti(d,1), d);
    2334       11395 :   return expIr(divri(mulri(mppi(prec), n), d));
    2335             : }
    2336             : /* exp(i Pi z), z a t_INT or t_FRAC */
    2337             : static GEN
    2338       39697 : expIPiQ(GEN z, long prec)
    2339             : {
    2340       39697 :   if (typ(z) == t_INT) return mpodd(z)? gen_m1: gen_1;
    2341        1988 :   return expIPifrac(z, prec);
    2342             : }
    2343             : 
    2344             : /* convert power of 2 t_REAL to rational */
    2345             : static GEN
    2346        9276 : real2nQ(GEN x)
    2347             : {
    2348        9276 :   long e = expo(x);
    2349             :   GEN z;
    2350        9276 :   if (e < 0)
    2351        2879 :     z = mkfrac(signe(x) < 0? gen_m1: gen_1, int2n(-e));
    2352             :   else
    2353             :   {
    2354        6397 :     z = int2n(e);
    2355        6397 :     if (signe(x) < 0) togglesign_safe(&z);
    2356             :   }
    2357        9276 :   return z;
    2358             : }
    2359             : /* x a real number */
    2360             : GEN
    2361      183756 : expIPiR(GEN x, long prec)
    2362             : {
    2363      183756 :   if (typ(x) == t_REAL && absrnz_equal2n(x)) x = real2nQ(x);
    2364      183756 :   switch(typ(x))
    2365             :   {
    2366        3231 :     case t_INT:  return mpodd(x)? gen_m1: gen_1;
    2367        1952 :     case t_FRAC: return expIPifrac(x, prec);
    2368             :   }
    2369      178573 :   return expIr(mulrr(mppi(prec), x));
    2370             : }
    2371             : /* z a t_COMPLEX */
    2372             : GEN
    2373      376996 : expIPiC(GEN z, long prec)
    2374             : {
    2375             :   GEN pi, r, x, y;
    2376      376996 :   if (typ(z) != t_COMPLEX) return expIPiR(z, prec);
    2377      194268 :   x = gel(z,1);
    2378      194268 :   y = gel(z,2); if (gequal0(y)) return expIPiR(x, prec);
    2379      193240 :   pi = mppi(prec);
    2380      193240 :   r = gmul(pi, y); togglesign(r); r = mpexp(r); /* exp(-pi y) */
    2381      193240 :   if (typ(x) == t_REAL && absrnz_equal2n(x)) x = real2nQ(x);
    2382      193240 :   switch(typ(x))
    2383             :   {
    2384       62330 :     case t_INT: if (mpodd(x)) togglesign(r);
    2385       62330 :                 return r;
    2386       11112 :     case t_FRAC: return gmul(r, expIPifrac(x, prec));
    2387             :   }
    2388      119798 :   return gmul(r, expIr(mulrr(pi, x)));
    2389             : }
    2390             : /* exp(I x y), more efficient for x in R, y pure imaginary */
    2391             : GEN
    2392      441701 : expIxy(GEN x, GEN y, long prec) { return gexp(gmul(x, mulcxI(y)), prec); }
    2393             : 
    2394             : static GEN
    2395       70949 : qq(GEN x, long prec)
    2396             : {
    2397       70949 :   long tx = typ(x);
    2398             :   GEN y;
    2399             : 
    2400       70949 :   if (is_scalar_t(tx))
    2401             :   {
    2402       70907 :     if (tx == t_PADIC) return x;
    2403       70893 :     x = upper_to_cx(x, &prec);
    2404       70879 :     return cxtoreal(expIPiC(gmul2n(x,1), prec)); /* e(x) */
    2405             :   }
    2406          42 :   if (! ( y = toser_i(x)) ) pari_err_TYPE("modular function", x);
    2407          42 :   return y;
    2408             : }
    2409             : 
    2410             : /* return (y * X^d) + x. Assume d > 0, x != 0, valser(x) = 0 */
    2411             : static GEN
    2412          21 : ser_addmulXn(GEN y, GEN x, long d)
    2413             : {
    2414          21 :   long i, lx, ly, l = valser(y) + d; /* > 0 */
    2415             :   GEN z;
    2416             : 
    2417          21 :   lx = lg(x);
    2418          21 :   ly = lg(y) + l; if (lx < ly) ly = lx;
    2419          21 :   if (l > lx-2) return gcopy(x);
    2420          21 :   z = cgetg(ly,t_SER);
    2421          77 :   for (i=2; i<=l+1; i++) gel(z,i) = gel(x,i);
    2422          70 :   for (   ; i < ly; i++) gel(z,i) = gadd(gel(x,i),gel(y,i-l));
    2423          21 :   z[1] = x[1]; return z;
    2424             : }
    2425             : 
    2426             : /* q a t_POL s.t. q(0) != 0, v > 0, Q = x^v*q; return \prod_i (1-Q^i) */
    2427             : static GEN
    2428          28 : RgXn_eta(GEN q, long v, long lim)
    2429             : {
    2430          28 :   pari_sp av = avma;
    2431             :   GEN qn, ps, y;
    2432             :   ulong vps, vqn, n;
    2433             : 
    2434          28 :   if (!degpol(q) && isint1(gel(q,2))) return eta_ZXn(v, lim+v);
    2435           7 :   y = qn = ps = pol_1(0);
    2436           7 :   vps = vqn = 0;
    2437           7 :   for(n = 0;; n++)
    2438           7 :   { /* qn = q^n,  ps = (-1)^n q^(n(3n+1)/2),
    2439             :      * vps, vqn valuation of ps, qn HERE */
    2440          14 :     pari_sp av2 = avma;
    2441          14 :     ulong vt = vps + 2*vqn + v; /* valuation of t at END of loop body */
    2442             :     long k1, k2;
    2443             :     GEN t;
    2444          14 :     vqn += v; vps = vt + vqn; /* valuation of qn, ps at END of body */
    2445          14 :     k1 = lim + v - vt + 1;
    2446          14 :     k2 = k1 - vqn; /* = lim + v - vps + 1 */
    2447          14 :     if (k1 <= 0) break;
    2448          14 :     t = RgX_mul(q, RgX_sqr(qn));
    2449          14 :     t = RgXn_red_shallow(t, k1);
    2450          14 :     t = RgX_mul(ps,t);
    2451          14 :     t = RgXn_red_shallow(t, k1);
    2452          14 :     t = RgX_neg(t); /* t = (-1)^(n+1) q^(n(3n+1)/2 + 2n+1) */
    2453          14 :     t = gerepileupto(av2, t);
    2454          14 :     y = RgX_addmulXn_shallow(t, y, vt);
    2455          14 :     if (k2 <= 0) break;
    2456             : 
    2457           7 :     qn = RgX_mul(qn,q);
    2458           7 :     ps = RgX_mul(t,qn);
    2459           7 :     ps = RgXn_red_shallow(ps, k2);
    2460           7 :     y = RgX_addmulXn_shallow(ps, y, vps);
    2461             : 
    2462           7 :     if (gc_needed(av,1))
    2463             :     {
    2464           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"eta, n = %ld", n);
    2465           0 :       gerepileall(av, 3, &y, &qn, &ps);
    2466             :     }
    2467             :   }
    2468           7 :   return y;
    2469             : }
    2470             : 
    2471             : static GEN
    2472       45155 : inteta(GEN q)
    2473             : {
    2474       45155 :   long tx = typ(q);
    2475             :   GEN ps, qn, y;
    2476             : 
    2477       45155 :   y = gen_1; qn = gen_1; ps = gen_1;
    2478       45155 :   if (tx==t_PADIC)
    2479             :   {
    2480          28 :     if (valp(q) <= 0) pari_err_DOMAIN("eta", "v_p(q)", "<=",gen_0,q);
    2481             :     for(;;)
    2482          56 :     {
    2483          77 :       GEN t = gneg_i(gmul(ps,gmul(q,gsqr(qn))));
    2484          77 :       y = gadd(y,t); qn = gmul(qn,q); ps = gmul(t,qn);
    2485          77 :       t = y;
    2486          77 :       y = gadd(y,ps); if (gequal(t,y)) return y;
    2487             :     }
    2488             :   }
    2489             : 
    2490       45127 :   if (tx == t_SER)
    2491             :   {
    2492             :     ulong vps, vqn;
    2493          42 :     long l = lg(q), v, n;
    2494             :     pari_sp av;
    2495             : 
    2496          42 :     v = valser(q); /* handle valuation separately to avoid overflow */
    2497          42 :     if (v <= 0) pari_err_DOMAIN("eta", "v_p(q)", "<=",gen_0,q);
    2498          35 :     y = ser2pol_i(q, l); /* t_SER inefficient when input has low degree */
    2499          35 :     n = degpol(y);
    2500          35 :     if (n <= (l>>2))
    2501             :     {
    2502          28 :       GEN z = RgXn_eta(y, v, l-2);
    2503          28 :       setvarn(z, varn(y)); return RgX_to_ser(z, l+v);
    2504             :     }
    2505           7 :     q = leafcopy(q); av = avma;
    2506           7 :     setvalser(q, 0);
    2507           7 :     y = scalarser(gen_1, varn(q), l+v);
    2508           7 :     vps = vqn = 0;
    2509           7 :     for(n = 0;; n++)
    2510           7 :     { /* qn = q^n,  ps = (-1)^n q^(n(3n+1)/2) */
    2511          14 :       ulong vt = vps + 2*vqn + v;
    2512             :       long k;
    2513             :       GEN t;
    2514          14 :       t = gneg_i(gmul(ps,gmul(q,gsqr(qn))));
    2515             :       /* t = (-1)^(n+1) q^(n(3n+1)/2 + 2n+1) */
    2516          14 :       y = ser_addmulXn(t, y, vt);
    2517          14 :       vqn += v; vps = vt + vqn;
    2518          14 :       k = l+v - vps; if (k <= 2) return y;
    2519             : 
    2520           7 :       qn = gmul(qn,q); ps = gmul(t,qn);
    2521           7 :       y = ser_addmulXn(ps, y, vps);
    2522           7 :       setlg(q, k);
    2523           7 :       setlg(qn, k);
    2524           7 :       setlg(ps, k);
    2525           7 :       if (gc_needed(av,3))
    2526             :       {
    2527           0 :         if(DEBUGMEM>1) pari_warn(warnmem,"eta");
    2528           0 :         gerepileall(av, 3, &y, &qn, &ps);
    2529             :       }
    2530             :     }
    2531             :   }
    2532             :   {
    2533       45085 :     long l = -prec2nbits(precision(q));
    2534       45085 :     pari_sp av = avma;
    2535             : 
    2536             :     for(;;)
    2537      128489 :     {
    2538      173574 :       GEN t = gneg_i(gmul(ps,gmul(q,gsqr(qn))));
    2539             :       /* qn = q^n
    2540             :        * ps = (-1)^n q^(n(3n+1)/2)
    2541             :        * t = (-1)^(n+1) q^(n(3n+1)/2 + 2n+1) */
    2542      173574 :       y = gadd(y,t); qn = gmul(qn,q); ps = gmul(t,qn);
    2543      173574 :       y = gadd(y,ps);
    2544      173574 :       if (gexpo(ps)-gexpo(y) < l) return y;
    2545      128489 :       if (gc_needed(av,3))
    2546             :       {
    2547           0 :         if(DEBUGMEM>1) pari_warn(warnmem,"eta");
    2548           0 :         gerepileall(av, 3, &y, &qn, &ps);
    2549             :       }
    2550             :     }
    2551             :   }
    2552             : }
    2553             : 
    2554             : GEN
    2555          77 : eta(GEN x, long prec)
    2556             : {
    2557          77 :   pari_sp av = avma;
    2558          77 :   GEN z = inteta( qq(x,prec) );
    2559          49 :   if (typ(z) == t_SER) return gerepilecopy(av, z);
    2560          14 :   return gerepileupto(av, z);
    2561             : }
    2562             : 
    2563             : /* s(h,k) = sum(n = 1, k-1, (n/k)*(frac(h*n/k) - 1/2))
    2564             :  * Knuth's algorithm. h integer, k integer > 0, (h,k) = 1 */
    2565             : GEN
    2566        6258 : sumdedekind_coprime(GEN h, GEN k)
    2567             : {
    2568        6258 :   pari_sp av = avma;
    2569             :   GEN s2, s1, p, pp;
    2570             :   long s;
    2571        6258 :   if (lgefint(k) == 3 && uel(k,2) <= (2*(ulong)LONG_MAX) / 3)
    2572             :   {
    2573        6251 :     ulong kk = k[2], hh = umodiu(h, kk);
    2574             :     long s1, s2;
    2575             :     GEN v;
    2576        6251 :     if (signe(k) < 0) { k = negi(k); hh = Fl_neg(hh, kk); }
    2577        6251 :     v = u_sumdedekind_coprime(hh, kk);
    2578        6251 :     s1 = v[1]; s2 = v[2];
    2579        6251 :     return gerepileupto(av, gdiv(addis(mulis(k,s1), s2), muluu(12, kk)));
    2580             :   }
    2581           7 :   s = 1;
    2582           7 :   s1 = gen_0; p = gen_1; pp = gen_0;
    2583           7 :   s2 = h = modii(h, k);
    2584          35 :   while (signe(h)) {
    2585          28 :     GEN r, nexth, a = dvmdii(k, h, &nexth);
    2586          28 :     if (is_pm1(h)) s2 = s == 1? addii(s2, p): subii(s2, p);
    2587          28 :     s1 = s == 1? addii(s1, a): subii(s1, a);
    2588          28 :     s = -s;
    2589          28 :     k = h; h = nexth;
    2590          28 :     r = addii(mulii(a,p), pp); pp = p; p = r;
    2591             :   }
    2592             :   /* at this point p = original k */
    2593           7 :   if (s == -1) s1 = subiu(s1, 3);
    2594           7 :   return gerepileupto(av, gdiv(addii(mulii(p,s1), s2), muliu(p,12)));
    2595             : }
    2596             : /* as above, for ulong arguments.
    2597             :  * k integer > 0, 0 <= h < k, (h,k) = 1. Returns [s1,s2] such that
    2598             :  * s(h,k) = (s2 + k s1) / (12k). Requires max(h + k/2, k) < LONG_MAX
    2599             :  * to avoid overflow, in particular k <= LONG_MAX * 2/3 is fine */
    2600             : GEN
    2601        6251 : u_sumdedekind_coprime(long h, long k)
    2602             : {
    2603        6251 :   long s = 1, s1 = 0, s2 = h, p = 1, pp = 0;
    2604       11424 :   while (h) {
    2605        5173 :     long r, nexth = k % h, a = k / h; /* a >= 1, a >= 2 if h = 1 */
    2606        5173 :     if (h == 1) s2 += p * s; /* occurs exactly once, last step */
    2607        5173 :     s1 += a * s;
    2608        5173 :     s = -s;
    2609        5173 :     k = h; h = nexth;
    2610        5173 :     r = a*p + pp; pp = p; p = r; /* p >= pp >= 0 */
    2611             :   }
    2612             :   /* in the above computation, p increases from 1 to original k,
    2613             :    * -k/2 <= s2 <= h + k/2, and |s1| <= k */
    2614        6251 :   if (s < 0) s1 -= 3; /* |s1| <= k+3 ? */
    2615             :   /* But in fact, |s2 + p s1| <= k^2 + 1/2 - 3k; if (s < 0), we have
    2616             :    * |s2| <= k/2 and it follows that |s1| < k here as well */
    2617             :   /* p = k; s(h,k) = (s2 + p s1)/12p. */
    2618        6251 :   return mkvecsmall2(s1, s2);
    2619             : }
    2620             : GEN
    2621          28 : sumdedekind(GEN h, GEN k)
    2622             : {
    2623          28 :   pari_sp av = avma;
    2624             :   GEN d;
    2625          28 :   if (typ(h) != t_INT) pari_err_TYPE("sumdedekind",h);
    2626          28 :   if (typ(k) != t_INT) pari_err_TYPE("sumdedekind",k);
    2627          28 :   d = gcdii(h,k);
    2628          28 :   if (!is_pm1(d))
    2629             :   {
    2630           7 :     h = diviiexact(h, d);
    2631           7 :     k = diviiexact(k, d);
    2632             :   }
    2633          28 :   return gerepileupto(av, sumdedekind_coprime(h,k));
    2634             : }
    2635             : 
    2636             : /* eta(x); assume Im x >> 0 (e.g. x in SL2's standard fundamental domain) */
    2637             : static GEN
    2638       45997 : eta_reduced(GEN x, long prec)
    2639             : {
    2640       45997 :   GEN z = expIPiC(gdivgu(x, 12), prec); /* e(x/24) */
    2641       45997 :   if (24 * gexpo(z) >= -prec2nbits(prec))
    2642       45064 :     z = gmul(z, inteta( gpowgs(z,24) ));
    2643       45997 :   return z;
    2644             : }
    2645             : 
    2646             : /* x = U.z (flag = 1), or x = U^(-1).z (flag = 0)
    2647             :  * Return [s,t] such that eta(z) = eta(x) * sqrt(s) * exp(I Pi t) */
    2648             : static GEN
    2649       46011 : eta_correction(GEN x, GEN U, long flag)
    2650             : {
    2651             :   GEN a,b,c,d, s,t;
    2652             :   long sc;
    2653       46011 :   a = gcoeff(U,1,1);
    2654       46011 :   b = gcoeff(U,1,2);
    2655       46011 :   c = gcoeff(U,2,1);
    2656       46011 :   d = gcoeff(U,2,2);
    2657             :   /* replace U by U^(-1) */
    2658       46011 :   if (flag) {
    2659       37639 :     swap(a,d);
    2660       37639 :     togglesign_safe(&b);
    2661       37639 :     togglesign_safe(&c);
    2662             :   }
    2663       46011 :   sc = signe(c);
    2664       46011 :   if (!sc) {
    2665       39781 :     if (signe(d) < 0) togglesign_safe(&b);
    2666       39781 :     s = gen_1;
    2667       39781 :     t = uutoQ(umodiu(b, 24), 12);
    2668             :   } else {
    2669        6230 :     if (sc < 0) {
    2670        1722 :       togglesign_safe(&a);
    2671        1722 :       togglesign_safe(&b);
    2672        1722 :       togglesign_safe(&c);
    2673        1722 :       togglesign_safe(&d);
    2674             :     } /* now c > 0 */
    2675        6230 :     s = mulcxmI(gadd(gmul(c,x), d));
    2676        6230 :     t = gadd(gdiv(addii(a,d),muliu(c,12)), sumdedekind_coprime(negi(d),c));
    2677             :     /* correction : exp(I Pi (((a+d)/12c) + s(-d,c)) ) sqrt(-i(cx+d))  */
    2678             :   }
    2679       46011 :   return mkvec2(s, t);
    2680             : }
    2681             : 
    2682             : /* returns the true value of eta(x) for Im(x) > 0, using reduction to
    2683             :  * standard fundamental domain */
    2684             : GEN
    2685       37569 : trueeta(GEN x, long prec)
    2686             : {
    2687       37569 :   pari_sp av = avma;
    2688             :   GEN U, st, s, t;
    2689             : 
    2690       37569 :   if (!is_scalar_t(typ(x))) pari_err_TYPE("trueeta",x);
    2691       37569 :   x = upper_to_cx(x, &prec);
    2692       37569 :   x = cxredsl2(x, &U);
    2693       37569 :   st = eta_correction(x, U, 1);
    2694       37569 :   x = eta_reduced(x, prec);
    2695       37569 :   s = gel(st, 1);
    2696       37569 :   t = gel(st, 2);
    2697       37569 :   x = gmul(x, expIPiQ(t, prec));
    2698       37569 :   if (s != gen_1) x = gmul(x, gsqrt(s, prec));
    2699       37569 :   return gerepileupto(av, x);
    2700             : }
    2701             : 
    2702             : GEN
    2703          77 : eta0(GEN x, long flag,long prec)
    2704          77 : { return flag? trueeta(x,prec): eta(x,prec); }
    2705             : 
    2706             : /* eta(q) = 1 + \sum_{n>0} (-1)^n * (q^(n(3n-1)/2) + q^(n(3n+1)/2)) */
    2707             : static GEN
    2708           7 : ser_eta(long prec)
    2709             : {
    2710           7 :   GEN e = cgetg(prec+2, t_SER), ed = e+2;
    2711             :   long n, j;
    2712           7 :   e[1] = evalsigne(1)|_evalvalser(0)|evalvarn(0);
    2713           7 :   gel(ed,0) = gen_1;
    2714         483 :   for (n = 1; n < prec; n++) gel(ed,n) = gen_0;
    2715          49 :   for (n = 1, j = 0; n < prec; n++)
    2716             :   {
    2717             :     GEN s;
    2718          49 :     j += 3*n-2; /* = n*(3*n-1) / 2 */;
    2719          49 :     if (j >= prec) break;
    2720          42 :     s = odd(n)? gen_m1: gen_1;
    2721          42 :     gel(ed, j) = s;
    2722          42 :     if (j+n >= prec) break;
    2723          42 :     gel(ed, j+n) = s;
    2724             :   }
    2725           7 :   return e;
    2726             : }
    2727             : 
    2728             : static GEN
    2729         476 : coeffEu(GEN fa)
    2730             : {
    2731         476 :   pari_sp av = avma;
    2732         476 :   return gerepileuptoint(av, mului(65520, usumdivk_fact(fa,11)));
    2733             : }
    2734             : /* E12 = 1 + q*E/691 */
    2735             : static GEN
    2736           7 : ser_E(long prec)
    2737             : {
    2738           7 :   GEN e = cgetg(prec+2, t_SER), ed = e+2;
    2739           7 :   GEN F = vecfactoru_i(2, prec); /* F[n] = factoru(n+1) */
    2740             :   long n;
    2741           7 :   e[1] = evalsigne(1)|_evalvalser(0)|evalvarn(0);
    2742           7 :   gel(ed,0) = utoipos(65520);
    2743         483 :   for (n = 1; n < prec; n++) gel(ed,n) = coeffEu(gel(F,n));
    2744           7 :   return e;
    2745             : }
    2746             : /* j = E12/Delta + 432000/691, E12 = 1 + q*E/691 */
    2747             : static GEN
    2748           7 : ser_j2(long prec, long v)
    2749             : {
    2750           7 :   pari_sp av = avma;
    2751           7 :   GEN iD = gpowgs(ginv(ser_eta(prec)), 24); /* q/Delta */
    2752           7 :   GEN J = gmul(ser_E(prec), iD);
    2753           7 :   setvalser(iD,-1); /* now 1/Delta */
    2754           7 :   J = gadd(gdivgu(J, 691), iD);
    2755           7 :   J = gerepileupto(av, J);
    2756           7 :   if (prec > 1) gel(J,3) = utoipos(744);
    2757           7 :   setvarn(J,v); return J;
    2758             : }
    2759             : 
    2760             : /* j(q) = \sum_{n >= -1} c(n)q^n,
    2761             :  * \sum_{n = -1}^{N-1} c(n) (-10n \sigma_3(N-n) + 21 \sigma_5(N-n))
    2762             :  * = c(N) (N+1)/24 */
    2763             : static GEN
    2764          14 : ser_j(long prec, long v)
    2765             : {
    2766             :   GEN j, J, S3, S5, F;
    2767             :   long i, n;
    2768          14 :   if (prec > 64) return ser_j2(prec, v);
    2769           7 :   S3 = cgetg(prec+1, t_VEC);
    2770           7 :   S5 = cgetg(prec+1,t_VEC);
    2771           7 :   F = vecfactoru_i(1, prec);
    2772          35 :   for (n = 1; n <= prec; n++)
    2773             :   {
    2774          28 :     GEN fa = gel(F,n);
    2775          28 :     gel(S3,n) = mului(10, usumdivk_fact(fa,3));
    2776          28 :     gel(S5,n) = mului(21, usumdivk_fact(fa,5));
    2777             :   }
    2778           7 :   J = cgetg(prec+2, t_SER),
    2779           7 :   J[1] = evalvarn(v)|evalsigne(1)|evalvalser(-1);
    2780           7 :   j = J+3;
    2781           7 :   gel(j,-1) = gen_1;
    2782           7 :   gel(j,0) = utoipos(744);
    2783           7 :   gel(j,1) = utoipos(196884);
    2784          21 :   for (n = 2; n < prec; n++)
    2785             :   {
    2786          14 :     pari_sp av = avma;
    2787          14 :     GEN c, s3 = gel(S3,n+1), s5 = gel(S5,n+1);
    2788          14 :     c = addii(s3, s5);
    2789          49 :     for (i = 0; i < n; i++)
    2790             :     {
    2791          35 :       s3 = gel(S3,n-i); s5 = gel(S5,n-i);
    2792          35 :       c = addii(c, mulii(gel(j,i), subii(s5, mului(i,s3))));
    2793             :     }
    2794          14 :     gel(j,n) = gerepileuptoint(av, diviuexact(muliu(c,24), n+1));
    2795             :   }
    2796           7 :   return J;
    2797             : }
    2798             : 
    2799             : GEN
    2800          42 : jell(GEN x, long prec)
    2801             : {
    2802          42 :   long tx = typ(x);
    2803          42 :   pari_sp av = avma;
    2804             :   GEN q, h, U;
    2805             : 
    2806          42 :   if (!is_scalar_t(tx))
    2807             :   {
    2808             :     long v;
    2809          21 :     if (gequalX(x)) return ser_j(precdl, varn(x));
    2810          21 :     q = toser_i(x); if (!q) pari_err_TYPE("ellj",x);
    2811          14 :     v = fetch_var_higher();
    2812          14 :     h = ser_j(lg(q)-2, v);
    2813          14 :     h = gsubst(h, v, q);
    2814          14 :     delete_var(); return gerepileupto(av, h);
    2815             :   }
    2816          21 :   if (tx == t_PADIC)
    2817             :   {
    2818           7 :     GEN p2, p1 = gdiv(inteta(gsqr(x)), inteta(x));
    2819           7 :     p1 = gmul2n(gsqr(p1),1);
    2820           7 :     p1 = gmul(x,gpowgs(p1,12));
    2821           7 :     p2 = gaddsg(768,gadd(gsqr(p1),gdivsg(4096,p1)));
    2822           7 :     p1 = gmulsg(48,p1);
    2823           7 :     return gerepileupto(av, gadd(p2,p1));
    2824             :   }
    2825             :   /* Let h = Delta(2x) / Delta(x), then j(x) = (1 + 256h)^3 / h */
    2826          14 :   x = upper_to_cx(x, &prec);
    2827           7 :   x = cxredsl2(x, &U); /* forget about Ua : j has weight 0 */
    2828             :   { /* cf eta_reduced, raised to power 24
    2829             :      * Compute
    2830             :      *   t = (inteta(q(2x)) / inteta(q(x))) ^ 24;
    2831             :      * then
    2832             :      *   h = t * (q(2x) / q(x) = t * q(x);
    2833             :      * but inteta(q) costly and useless if expo(q) << 1  => inteta(q) = 1.
    2834             :      * log_2 ( exp(-2Pi Im tau) ) < -prec2nbits(prec)
    2835             :      * <=> Im tau > prec2nbits(prec) * log(2) / 2Pi */
    2836           7 :     long C = (long)prec2nbits_mul(prec, M_LN2/(2*M_PI));
    2837           7 :     q = expIPiC(gmul2n(x,1), prec); /* e(x) */
    2838           7 :     if (gcmpgs(gel(x,2), C) > 0) /* eta(q(x)) = 1 : no need to compute q(2x) */
    2839           0 :       h = q;
    2840             :     else
    2841             :     {
    2842           7 :       GEN t = gdiv(inteta(gsqr(q)), inteta(q));
    2843           7 :       h = gmul(q, gpowgs(t, 24));
    2844             :     }
    2845             :   }
    2846             :   /* real_1 important ! gaddgs(, 1) could increase the accuracy ! */
    2847           7 :   return gerepileupto(av, gdiv(gpowgs(gadd(gmul2n(h,8), real_1(prec)), 3), h));
    2848             : }
    2849             : 
    2850             : static GEN
    2851        8372 : to_form(GEN a, GEN w, GEN C, GEN D)
    2852        8372 : { return mkqfb(a, w, diviiexact(C, a), D); }
    2853             : static GEN
    2854        8372 : form_to_quad(GEN f, GEN sqrtD)
    2855             : {
    2856        8372 :   long a = itos(gel(f,1)), a2 = a << 1;
    2857        8372 :   GEN b = gel(f,2);
    2858        8372 :   return mkcomplex(gdivgs(b, -a2), gdivgs(sqrtD, a2));
    2859             : }
    2860             : static GEN
    2861        8372 : eta_form(GEN f, GEN sqrtD, GEN *s_t, long prec)
    2862             : {
    2863        8372 :   GEN U, t = form_to_quad(redimagsl2(f, &U), sqrtD);
    2864        8372 :   *s_t = eta_correction(t, U, 0);
    2865        8372 :   return eta_reduced(t, prec);
    2866             : }
    2867             : 
    2868             : /* eta(t/p)eta(t/q) / (eta(t)eta(t/pq)), t = (-w + sqrt(D)) / 2a */
    2869             : GEN
    2870        2093 : double_eta_quotient(GEN a, GEN w, GEN D, long p, long q, GEN pq, GEN sqrtD)
    2871             : {
    2872        2093 :   GEN C = shifti(subii(sqri(w), D), -2);
    2873             :   GEN d, t, z, zp, zq, zpq, s_t, s_tp, s_tpq, s, sp, spq;
    2874        2093 :   long prec = realprec(sqrtD);
    2875             : 
    2876        2093 :   z = eta_form(to_form(a, w, C, D), sqrtD, &s_t, prec);
    2877        2093 :   s = gel(s_t, 1);
    2878        2093 :   zp = eta_form(to_form(mului(p, a), w, C, D), sqrtD, &s_tp, prec);
    2879        2093 :   sp = gel(s_tp, 1);
    2880        2093 :   zpq = eta_form(to_form(mulii(pq, a), w, C, D), sqrtD, &s_tpq, prec);
    2881        2093 :   spq = gel(s_tpq, 1);
    2882        2093 :   if (p == q) {
    2883           0 :     z = gdiv(gsqr(zp), gmul(z, zpq));
    2884           0 :     t = gsub(gmul2n(gel(s_tp,2), 1),
    2885           0 :              gadd(gel(s_t,2), gel(s_tpq,2)));
    2886           0 :     if (sp != gen_1) z = gmul(z, sp);
    2887             :   } else {
    2888             :     GEN s_tq, sq;
    2889        2093 :     zq = eta_form(to_form(mului(q, a), w, C, D), sqrtD, &s_tq, prec);
    2890        2093 :     sq = gel(s_tq, 1);
    2891        2093 :     z = gdiv(gmul(zp, zq), gmul(z, zpq));
    2892        2093 :     t = gsub(gadd(gel(s_tp,2), gel(s_tq,2)),
    2893        2093 :              gadd(gel(s_t,2), gel(s_tpq,2)));
    2894        2093 :     if (sp != gen_1) z = gmul(z, gsqrt(sp, prec));
    2895        2093 :     if (sq != gen_1) z = gmul(z, gsqrt(sq, prec));
    2896             :   }
    2897        2093 :   d = NULL;
    2898        2093 :   if (s != gen_1) d = gsqrt(s, prec);
    2899        2093 :   if (spq != gen_1) {
    2900        2065 :     GEN x = gsqrt(spq, prec);
    2901        2065 :     d = d? gmul(d, x): x;
    2902             :   }
    2903        2093 :   if (d) z = gdiv(z, d);
    2904        2093 :   return gmul(z, expIPiQ(t, prec));
    2905             : }
    2906             : 
    2907             : typedef struct { GEN u; long v, t; } cxanalyze_t;
    2908             : 
    2909             : /* Check whether a t_COMPLEX, t_REAL or t_INT z != 0 can be written as
    2910             :  * z = u * 2^(v/2) * exp(I Pi/4 t), u > 0, v = 0,1 and -3 <= t <= 4.
    2911             :  * Allow z t_INT/t_REAL to simplify handling of eta_correction() output */
    2912             : static int
    2913          70 : cxanalyze(cxanalyze_t *T, GEN z)
    2914             : {
    2915             :   GEN a, b;
    2916             :   long ta, tb;
    2917             : 
    2918          70 :   T->v = 0;
    2919          70 :   if (is_intreal_t(typ(z)))
    2920             :   {
    2921          63 :     T->u = mpabs_shallow(z);
    2922          63 :     T->t = signe(z) < 0? 4: 0;
    2923          63 :     return 1;
    2924             :   }
    2925           7 :   a = gel(z,1); ta = typ(a);
    2926           7 :   b = gel(z,2); tb = typ(b);
    2927             : 
    2928           7 :   T->t = 0;
    2929           7 :   if (ta == t_INT && !signe(a))
    2930             :   {
    2931           0 :     T->u = R_abs_shallow(b);
    2932           0 :     T->t = gsigne(b) < 0? -2: 2;
    2933           0 :     return 1;
    2934             :   }
    2935           7 :   if (tb == t_INT && !signe(b))
    2936             :   {
    2937           0 :     T->u = R_abs_shallow(a);
    2938           0 :     T->t = gsigne(a) < 0? 4: 0;
    2939           0 :     return 1;
    2940             :   }
    2941           7 :   if (ta != tb || ta == t_REAL) { T->u = z; return 0; }
    2942             :   /* a,b both non zero, both t_INT or t_FRAC */
    2943           7 :   if (ta == t_INT)
    2944             :   {
    2945           7 :     if (!absequalii(a, b)) return 0;
    2946           7 :     T->u = absi_shallow(a);
    2947           7 :     T->v = 1;
    2948           7 :     if (signe(a) == signe(b))
    2949           0 :     { T->t = signe(a) < 0? -3: 1; }
    2950             :     else
    2951           7 :     { T->t = signe(a) < 0? 3: -1; }
    2952             :   }
    2953             :   else
    2954             :   {
    2955           0 :     if (!absequalii(gel(a,2), gel(b,2)) || !absequalii(gel(a,1),gel(b,1)))
    2956           0 :       return 0;
    2957           0 :     T->u = absfrac_shallow(a);
    2958           0 :     T->v = 1;
    2959           0 :     a = gel(a,1);
    2960           0 :     b = gel(b,1);
    2961           0 :     if (signe(a) == signe(b))
    2962           0 :     { T->t = signe(a) < 0? -3: 1; }
    2963             :     else
    2964           0 :     { T->t = signe(a) < 0? 3: -1; }
    2965             :   }
    2966           7 :   return 1;
    2967             : }
    2968             : 
    2969             : /* z * sqrt(st_b) / sqrt(st_a) exp(I Pi (t + t0)). Assume that
    2970             :  * sqrt2 = gsqrt(gen_2, prec) or NULL */
    2971             : static GEN
    2972          35 : apply_eta_correction(GEN z, GEN st_a, GEN st_b, GEN t0, GEN sqrt2, long prec)
    2973             : {
    2974          35 :   GEN t, s_a = gel(st_a, 1), s_b = gel(st_b, 1);
    2975             :   cxanalyze_t Ta, Tb;
    2976             :   int ca, cb;
    2977             : 
    2978          35 :   t = gsub(gel(st_b,2), gel(st_a,2));
    2979          35 :   if (t0 != gen_0) t = gadd(t, t0);
    2980          35 :   ca = cxanalyze(&Ta, s_a);
    2981          35 :   cb = cxanalyze(&Tb, s_b);
    2982          35 :   if (ca || cb)
    2983          35 :   { /* compute sqrt(s_b) / sqrt(s_a) in a more efficient way:
    2984             :      * sb = ub sqrt(2)^vb exp(i Pi/4 tb) */
    2985          35 :     GEN u = gdiv(Tb.u,Ta.u);
    2986          35 :     switch(Tb.v - Ta.v)
    2987             :     {
    2988           0 :       case -1: u = gmul2n(u,-1); /* fall through: write 1/sqrt2 = sqrt2/2 */
    2989           7 :       case 1: u = gmul(u, sqrt2? sqrt2: sqrtr_abs(real2n(1, prec)));
    2990             :     }
    2991          35 :     if (!isint1(u)) z = gmul(z, gsqrt(u, prec));
    2992          35 :     t = gadd(t, gmul2n(stoi(Tb.t - Ta.t), -3));
    2993             :   }
    2994             :   else
    2995             :   {
    2996           0 :     z = gmul(z, gsqrt(s_b, prec));
    2997           0 :     z = gdiv(z, gsqrt(s_a, prec));
    2998             :   }
    2999          35 :   return gmul(z, expIPiQ(t, prec));
    3000             : }
    3001             : 
    3002             : /* sqrt(2) eta(2x) / eta(x) */
    3003             : GEN
    3004           7 : weberf2(GEN x, long prec)
    3005             : {
    3006           7 :   pari_sp av = avma;
    3007             :   GEN z, sqrt2, a,b, Ua,Ub, st_a,st_b;
    3008             : 
    3009           7 :   x = upper_to_cx(x, &prec);
    3010           7 :   a = cxredsl2(x, &Ua);
    3011           7 :   b = cxredsl2(gmul2n(x,1), &Ub);
    3012           7 :   if (gequal(a,b)) /* not infrequent */
    3013           0 :     z = gen_1;
    3014             :   else
    3015           7 :     z = gdiv(eta_reduced(b,prec), eta_reduced(a,prec));
    3016           7 :   st_a = eta_correction(a, Ua, 1);
    3017           7 :   st_b = eta_correction(b, Ub, 1);
    3018           7 :   sqrt2 = sqrtr_abs(real2n(1, prec));
    3019           7 :   z = apply_eta_correction(z, st_a, st_b, gen_0, sqrt2, prec);
    3020           7 :   return gerepileupto(av, gmul(z, sqrt2));
    3021             : }
    3022             : 
    3023             : /* eta(x/2) / eta(x) */
    3024             : GEN
    3025          14 : weberf1(GEN x, long prec)
    3026             : {
    3027          14 :   pari_sp av = avma;
    3028             :   GEN z, a,b, Ua,Ub, st_a,st_b;
    3029             : 
    3030          14 :   x = upper_to_cx(x, &prec);
    3031          14 :   a = cxredsl2(x, &Ua);
    3032          14 :   b = cxredsl2(gmul2n(x,-1), &Ub);
    3033          14 :   if (gequal(a,b)) /* not infrequent */
    3034           0 :     z = gen_1;
    3035             :   else
    3036          14 :     z = gdiv(eta_reduced(b,prec), eta_reduced(a,prec));
    3037          14 :   st_a = eta_correction(a, Ua, 1);
    3038          14 :   st_b = eta_correction(b, Ub, 1);
    3039          14 :   z = apply_eta_correction(z, st_a, st_b, gen_0, NULL, prec);
    3040          14 :   return gerepileupto(av, z);
    3041             : }
    3042             : /* exp(-I*Pi/24) * eta((x+1)/2) / eta(x) */
    3043             : GEN
    3044          14 : weberf(GEN x, long prec)
    3045             : {
    3046          14 :   pari_sp av = avma;
    3047             :   GEN z, t0, a,b, Ua,Ub, st_a,st_b;
    3048          14 :   x = upper_to_cx(x, &prec);
    3049          14 :   a = cxredsl2(x, &Ua);
    3050          14 :   b = cxredsl2(gmul2n(gaddgs(x,1),-1), &Ub);
    3051          14 :   if (gequal(a,b)) /* not infrequent */
    3052           7 :     z = gen_1;
    3053             :   else
    3054           7 :     z = gdiv(eta_reduced(b,prec), eta_reduced(a,prec));
    3055          14 :   st_a = eta_correction(a, Ua, 1);
    3056          14 :   st_b = eta_correction(b, Ub, 1);
    3057          14 :   t0 = mkfrac(gen_m1, utoipos(24));
    3058          14 :   z = apply_eta_correction(z, st_a, st_b, t0, NULL, prec);
    3059          14 :   if (typ(z) == t_COMPLEX && isexactzero(real_i(x)))
    3060           0 :     z = gerepilecopy(av, gel(z,1));
    3061             :   else
    3062          14 :     z = gerepileupto(av, z);
    3063          14 :   return z;
    3064             : }
    3065             : GEN
    3066          35 : weber0(GEN x, long flag,long prec)
    3067             : {
    3068          35 :   switch(flag)
    3069             :   {
    3070          14 :     case 0: return weberf(x,prec);
    3071          14 :     case 1: return weberf1(x,prec);
    3072           7 :     case 2: return weberf2(x,prec);
    3073           0 :     default: pari_err_FLAG("weber");
    3074             :   }
    3075             :   return NULL; /* LCOV_EXCL_LINE */
    3076             : }
    3077             : 
    3078             : /* check |q| < 1 */
    3079             : static GEN
    3080          21 : check_unit_disc(const char *fun, GEN q, long prec)
    3081             : {
    3082          21 :   GEN Q = gtofp(q, prec), Qlow;
    3083          21 :   Qlow = (prec > LOWDEFAULTPREC)? gtofp(Q,LOWDEFAULTPREC): Q;
    3084          21 :   if (gcmp(gnorm(Qlow), gen_1) >= 0)
    3085           0 :     pari_err_DOMAIN(fun, "abs(q)", ">=", gen_1, q);
    3086          21 :   return Q;
    3087             : }
    3088             : 
    3089             : GEN
    3090          14 : theta(GEN q, GEN z, long prec)
    3091             : {
    3092             :   long l, n;
    3093          14 :   pari_sp av = avma, av2;
    3094             :   GEN s, c, snz, cnz, s2z, c2z, ps, qn, y, zy, ps2, k, zold;
    3095             : 
    3096          14 :   l = precision(q);
    3097          14 :   n = precision(z); if (n && n < l) l = n;
    3098          14 :   if (l) prec = l;
    3099          14 :   z = gtofp(z, prec);
    3100          14 :   q = check_unit_disc("theta", q, prec);
    3101          14 :   zold = NULL; /* gcc -Wall */
    3102          14 :   zy = imag_i(z);
    3103          14 :   if (gequal0(zy)) k = gen_0;
    3104             :   else
    3105             :   {
    3106           7 :     GEN lq = glog(q,prec); k = roundr(divrr(zy, real_i(lq)));
    3107           7 :     if (signe(k)) { zold = z; z = gadd(z, mulcxmI(gmul(lq,k))); }
    3108             :   }
    3109          14 :   qn = gen_1;
    3110          14 :   ps2 = gsqr(q);
    3111          14 :   ps = gneg_i(ps2);
    3112          14 :   gsincos(z, &s, &c, prec);
    3113          14 :   s2z = gmul2n(gmul(s,c),1); /* sin 2z */
    3114          14 :   c2z = gsubgs(gmul2n(gsqr(c),1), 1); /* cos 2z */
    3115          14 :   snz = s;
    3116          14 :   cnz = c; y = s;
    3117          14 :   av2 = avma;
    3118          14 :   for (n = 3;; n += 2)
    3119         259 :   {
    3120             :     long e;
    3121         273 :     s = gadd(gmul(snz, c2z), gmul(cnz,s2z));
    3122         273 :     qn = gmul(qn,ps);
    3123         273 :     y = gadd(y, gmul(qn, s));
    3124         273 :     e = gexpo(s); if (e < 0) e = 0;
    3125         273 :     if (gexpo(qn) + e < -prec2nbits(prec)) break;
    3126             : 
    3127         259 :     ps = gmul(ps,ps2);
    3128         259 :     c = gsub(gmul(cnz, c2z), gmul(snz,s2z));
    3129         259 :     snz = s; /* sin nz */
    3130         259 :     cnz = c; /* cos nz */
    3131         259 :     if (gc_needed(av,2))
    3132             :     {
    3133           0 :       if (DEBUGMEM>1) pari_warn(warnmem,"theta (n = %ld)", n);
    3134           0 :       gerepileall(av2, 5, &snz, &cnz, &ps, &qn, &y);
    3135             :     }
    3136             :   }
    3137          14 :   if (signe(k))
    3138             :   {
    3139           7 :     y = gmul(y, gmul(powgi(q,sqri(k)),
    3140             :                      gexp(gmul(mulcxI(zold),shifti(k,1)), prec)));
    3141           7 :     if (mod2(k)) y = gneg_i(y);
    3142             :   }
    3143          14 :   return gerepileupto(av, gmul(y, gmul2n(gsqrt(gsqrt(q,prec),prec),1)));
    3144             : }
    3145             : 
    3146             : GEN
    3147           7 : thetanullk(GEN q, long k, long prec)
    3148             : {
    3149             :   long l, n;
    3150           7 :   pari_sp av = avma;
    3151             :   GEN p1, ps, qn, y, ps2;
    3152             : 
    3153           7 :   if (k < 0)
    3154           0 :     pari_err_DOMAIN("thetanullk", "k", "<", gen_0, stoi(k));
    3155           7 :   l = precision(q);
    3156           7 :   if (l) prec = l;
    3157           7 :   q = check_unit_disc("thetanullk", q, prec);
    3158             : 
    3159           7 :   if (!odd(k)) { set_avma(av); return gen_0; }
    3160           7 :   qn = gen_1;
    3161           7 :   ps2 = gsqr(q);
    3162           7 :   ps = gneg_i(ps2);
    3163           7 :   y = gen_1;
    3164           7 :   for (n = 3;; n += 2)
    3165         280 :   {
    3166             :     GEN t;
    3167         287 :     qn = gmul(qn,ps);
    3168         287 :     ps = gmul(ps,ps2);
    3169         287 :     t = gmul(qn, powuu(n, k)); y = gadd(y, t);
    3170         287 :     if (gexpo(t) < -prec2nbits(prec)) break;
    3171             :   }
    3172           7 :   p1 = gmul2n(gsqrt(gsqrt(q,prec),prec),1);
    3173           7 :   if (k&2) y = gneg_i(y);
    3174           7 :   return gerepileupto(av, gmul(p1, y));
    3175             : }
    3176             : 
    3177             : /* q2 = q^2 */
    3178             : static GEN
    3179       70872 : vecthetanullk_loop(GEN q2, long k, long prec)
    3180             : {
    3181       70872 :   GEN ps, qn = gen_1, y = const_vec(k, gen_1);
    3182       70872 :   pari_sp av = avma;
    3183       70872 :   const long bit = prec2nbits(prec);
    3184             :   long i, n;
    3185             : 
    3186       70872 :   if (gexpo(q2) < -2*bit) return y;
    3187       70872 :   ps = gneg_i(q2);
    3188       70872 :   for (n = 3;; n += 2)
    3189      359719 :   {
    3190      430591 :     GEN t = NULL/*-Wall*/, P = utoipos(n), N2 = sqru(n);
    3191      430591 :     qn = gmul(qn,ps);
    3192      430591 :     ps = gmul(ps,q2);
    3193     1291773 :     for (i = 1; i <= k; i++)
    3194             :     {
    3195      861182 :       t = gmul(qn, P); gel(y,i) = gadd(gel(y,i), t);
    3196      861182 :       P = mulii(P, N2);
    3197             :     }
    3198      430591 :     if (gexpo(t) < -bit) return y;
    3199      359719 :     if (gc_needed(av,2))
    3200             :     {
    3201           0 :       if (DEBUGMEM>1) pari_warn(warnmem,"vecthetanullk_loop, n = %ld",n);
    3202           0 :       gerepileall(av, 3, &qn, &ps, &y);
    3203             :     }
    3204             :   }
    3205             : }
    3206             : /* [d^i theta/dz^i(q, 0), i = 1, 3, .., 2*k - 1] */
    3207             : GEN
    3208           0 : vecthetanullk(GEN q, long k, long prec)
    3209             : {
    3210           0 :   long i, l = precision(q);
    3211           0 :   pari_sp av = avma;
    3212             :   GEN p1, y;
    3213             : 
    3214           0 :   if (l) prec = l;
    3215           0 :   q = check_unit_disc("vecthetanullk", q, prec);
    3216           0 :   y = vecthetanullk_loop(gsqr(q), k, prec);
    3217           0 :   p1 = gmul2n(gsqrt(gsqrt(q,prec),prec),1);
    3218           0 :   for (i = 2; i <= k; i += 2) gel(y,i) = gneg_i(gel(y,i));
    3219           0 :   return gerepileupto(av, gmul(p1, y));
    3220             : }
    3221             : 
    3222             : /* [d^i theta/dz^i(q, 0), i = 1, 3, .., 2*k - 1], q = exp(2iPi tau) */
    3223             : GEN
    3224           0 : vecthetanullk_tau(GEN tau, long k, long prec)
    3225             : {
    3226           0 :   long i, l = precision(tau);
    3227           0 :   pari_sp av = avma;
    3228             :   GEN q4, y;
    3229             : 
    3230           0 :   if (l) prec = l;
    3231           0 :   if (typ(tau) != t_COMPLEX || gsigne(gel(tau,2)) <= 0)
    3232           0 :     pari_err_DOMAIN("vecthetanullk_tau", "imag(tau)", "<=", gen_0, tau);
    3233           0 :   q4 = expIPiC(gmul2n(tau,-1), prec); /* q^(1/4) */
    3234           0 :   y = vecthetanullk_loop(gpowgs(q4,8), k, prec);
    3235           0 :   for (i = 2; i <= k; i += 2) gel(y,i) = gneg_i(gel(y,i));
    3236           0 :   return gerepileupto(av, gmul(gmul2n(q4,1), y));
    3237             : }
    3238             : 
    3239             : /* Return E_k(tau). Slow if tau is not in standard fundamental domain */
    3240             : GEN
    3241       71141 : cxEk(GEN tau, long k, long prec)
    3242             : {
    3243             :   pari_sp av;
    3244             :   GEN q, y, qn;
    3245       71141 :   long n, b, l = precision(tau);
    3246             : 
    3247       71141 :   if (l) prec = l;
    3248       71141 :   b = prec2nbits(prec);
    3249             :   /* sum n^(k-1) x^n <= x(1 + (k!-1)x) / (1-x)^k (cf Eulerian polynomials)
    3250             :    * S = \sum_{n > 0} n^(k-1) |q^n/(1-q^n)| <= x(1+(k!-1)x) / (1-x)^(k+1),
    3251             :    * where x = |q| = exp(-2Pi Im(tau)) < 1. Neglegt 2/zeta(1-k) * S if
    3252             :    * (2Pi)^k/(k-1)! x < 2^(-b-1) and k! x < 1. Use log2((2Pi)^k/(k-1)!) < 10 */
    3253       71141 :   if (gcmpgs(imag_i(tau), (M_LN2 / (2*M_PI)) * (b+1+10)) > 0)
    3254          17 :     return real_1(prec);
    3255       71124 :   if (k == 2)
    3256             :   { /* -theta^(3)(tau/2) / theta^(1)(tau/2). Assume that Im tau > 0 */
    3257       70872 :     y = vecthetanullk_loop(qq(tau,prec), 2, prec);
    3258       70872 :     return gdiv(gel(y,2), gel(y,1));
    3259             :   }
    3260         252 :   q = cxtoreal(expIPiC(gneg(gmul2n(tau,1)), prec));
    3261         252 :   av = avma; y = gen_0; qn = q;
    3262         252 :   for(n = 1;; n++)
    3263        2623 :   { /* compute y := sum_{n>0} n^(k-1) / (q^n-1) */
    3264        2875 :     GEN t = gdiv(powuu(n,k-1), gsubgs(qn,1));
    3265        2875 :     if (gequal0(t) || gexpo(t) <= - prec2nbits(prec) - 5) break;
    3266        2623 :     y = gadd(y, t);
    3267        2623 :     qn = gmul(q,qn);
    3268        2623 :     if (gc_needed(av,2))
    3269             :     {
    3270           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"elleisnum");
    3271           0 :       gerepileall(av, 2, &y,&qn);
    3272             :     }
    3273             :   }
    3274         252 :   return gadd(gen_1, gmul(y, gdiv(gen_2, szeta(1-k, prec))));
    3275             : }

Generated by: LCOV version 1.14