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.14.0 lcov report (development 27775-aca467eab2) Lines: 1684 1807 93.2 %
Date: 2022-07-03 07:33:15 Functions: 124 127 97.6 %
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      488848 : _abs(GEN x)
      36      488848 : { return gabs(gtofp(x,LOWDEFAULTPREC), LOWDEFAULTPREC); }
      37             : /* can we use asymptotic expansion ? */
      38             : static int
      39      405158 : bessel_asymp(GEN z, long bit)
      40      405158 : { return gcmpgs(_abs(z), (bit+10)/2) >= 0; }
      41             : 
      42             : /* Region I: 0 < Arg z <= Pi, II: -Pi < Arg z <= 0 */
      43             : static int
      44         532 : regI(GEN z)
      45             : {
      46         532 :   long s = gsigne(imag_i(z));
      47         532 :   return (s > 0 || (s == 0 && gsigne(real_i(z)) < 0)) ? 1 : 2;
      48             : }
      49             : /* Region 1: Re(z) >= 0, 2: Re(z) < 0, Im(z) >= 0, 3: Re(z) < 0, Im(z) < 0 */
      50             : static int
      51       82696 : regJ(GEN z)
      52             : {
      53       82696 :   if (gsigne(real_i(z)) >= 0) return 1;
      54         336 :   return gsigne(imag_i(z)) >= 0 ? 2 : 3;
      55             : }
      56             : 
      57             : /* Hankel's expansions:
      58             :  * a_k(n) = \prod_{0 <= j < k} (4n^2 - (2j+1)^2)
      59             :  * C(k)[n,z] = a_k(n) / (k! (8 z)^k)
      60             :  * A(z)  = exp(-z) sum_{k >= 0} C(k)
      61             :  * A(-z) = exp(z) sum_{k >= 0} (-1)^k C(k)
      62             :  * J_n(z) ~ [1] (A(z/i) / r + A(-z/i) r) / sqrt(2Pi z)
      63             :  *          [2] (A(z/i) r^3 + A(-z/i) r) / sqrt(2Pi z)
      64             :  *          [3] (A(z/i) / r + A(-z/i) / r^3) / sqrt(2Pi z)
      65             :  * Y_n(z) ~ [1] i(A(z/i) / r + A(-z/i) r) / sqrt(2Pi z)
      66             :  *          [2] i(A(z/i) (r^3-2/r) + A(-z/i) r) / sqrt(2Pi z)
      67             :  *          [3] i(-A(z/i)/r + A(-z/i)(2r-1/r^3)) / sqrt(2Pi z)
      68             :  * K_n(z) ~ A(z) Pi / sqrt(2 Pi z)
      69             :  * I_n(z) ~ [I] (A(-z) + r^2 A(z)) / sqrt(2 Pi z)
      70             :  *          [II](A(-z) + r^(-2) A(z)) / sqrt(2 Pi z) */
      71             : 
      72             : /* set [A(z), A(-z), exp((2*nu+1)*I*Pi/4)] */
      73             : static void
      74       83690 : hankel_ABr(GEN *pA, GEN *pB, GEN *pr, GEN n, GEN z, long bit)
      75             : {
      76       83690 :   GEN E, P, C, Q = gen_0, zi = ginv(gmul2n(z, 3));
      77       83690 :   GEN K = gaddgs(_abs(n), 1), n2 = gmul2n(gsqr(n),2);
      78       83690 :   long prec = nbits2prec(bit), B = bit + 4, m;
      79             : 
      80       83690 :   P = C = real_1_bit(bit);
      81       83690 :   for (m = 1;; m += 2)
      82             :   {
      83     5585602 :     C = gmul(C, gdivgu(gmul(gsub(n2, sqru(2*m - 1)), zi), m));
      84     5585602 :     Q = gadd(Q, C);
      85     5585602 :     C = gmul(C, gdivgu(gmul(gsub(n2, sqru(2*m + 1)), zi), m + 1));
      86     5585602 :     P = gadd(P, C);
      87     5585602 :     if (gexpo(C) < -B && gcmpgs(K, m) <= 0) break;
      88             :   }
      89       83690 :   E = gexp(z, prec);
      90       83690 :   *pA = gdiv(gadd(P, Q), E);
      91       83690 :   *pB = gmul(gsub(P, Q), E);
      92       83690 :   *pr = gexp(mulcxI(gmul(gaddgs(gmul2n(n,1), 1), Pi2n(-2, prec))), prec);
      93       83690 : }
      94             : 
      95             : /* sqrt(2*Pi*z) */
      96             : static GEN
      97       83690 : sqz(GEN z, long bit)
      98             : {
      99       83690 :   long prec = nbits2prec(bit);
     100       83690 :   return gsqrt(gmul(Pi2n(1, prec), z), prec);
     101             : }
     102             : 
     103             : static GEN
     104         462 : besskasymp(GEN nu, GEN z, long bit)
     105             : {
     106             :   GEN A, B, r;
     107         462 :   long prec = nbits2prec(bit);
     108         462 :   hankel_ABr(&A,&B,&r, nu, z, bit);
     109         462 :   return gdiv(gmul(A, mppi(prec)), sqz(z, bit));
     110             : }
     111             : 
     112             : static GEN
     113         532 : bessiasymp(GEN nu, GEN z, long bit)
     114             : {
     115             :   GEN A, B, r, R, r2;
     116         532 :   hankel_ABr(&A,&B,&r, nu, z, bit);
     117         532 :   r2 = gsqr(r);
     118         532 :   R = regI(z) == 1 ? gmul(A, r2) : gdiv(A, r2);
     119         532 :   return gdiv(gadd(B, R), sqz(z, bit));
     120             : }
     121             : 
     122             : static GEN
     123       82226 : bessjasymp(GEN nu, GEN z, long bit)
     124             : {
     125             :   GEN A, B, r, R;
     126       82226 :   long reg = regJ(z);
     127       82226 :   hankel_ABr(&A,&B,&r, nu, mulcxmI(z), bit);
     128       82226 :   if (reg == 1) R = gadd(gdiv(A, r), gmul(B, r));
     129         168 :   else if (reg == 2) R = gadd(gmul(A, gpowgs(r, 3)), gmul(B, r));
     130          56 :   else R = gadd(gdiv(A, r), gdiv(B, gpowgs(r, 3)));
     131       82226 :   return gdiv(R, sqz(z, bit));
     132             : }
     133             : 
     134             : static GEN
     135         470 : bessyasymp(GEN nu, GEN z, long bit)
     136             : {
     137             :   GEN A, B, r, R;
     138         470 :   long reg = regJ(z);
     139         470 :   hankel_ABr(&A,&B,&r, nu, mulcxmI(z), bit);
     140         470 :   if (reg == 1) R = gsub(gmul(B, r), gdiv(A, r));
     141         168 :   else if (reg == 2)
     142         112 :     R = gadd(gmul(A, gsub(gpowgs(r, 3), gmul2n(ginv(r), 1))), gmul(B, r));
     143             :   else
     144          56 :     R = gsub(gmul(B, gsub(gmul2n(r, 1), ginv(gpowgs(r, 3)))), gdiv(A, r));
     145         470 :   return gdiv(mulcxI(R), sqz(z, bit));
     146             : }
     147             : 
     148             : /* n! sum_{0 <= k <= m} x^k / (k!*(k+n)!) */
     149             : static GEN
     150      318172 : _jbessel(GEN n, GEN x, long m)
     151             : {
     152      318172 :   pari_sp av = avma;
     153      318172 :   GEN s = gen_1;
     154             :   long k;
     155             : 
     156   118000260 :   for (k = m; k >= 1; k--)
     157             :   {
     158   117682088 :     s = gaddsg(1, gdiv(gmul(x,s), gmulgu(gaddgs(n, k), k)));
     159   117682088 :     if (gc_needed(av,1))
     160             :     {
     161           0 :       if (DEBUGMEM>1) pari_warn(warnmem,"besselj");
     162           0 :       s = gerepileupto(av, s);
     163             :     }
     164             :   }
     165      318172 :   return s;
     166             : }
     167             : 
     168             : /* max(2, L * approximate solution to x log x = B) */
     169             : static long
     170      320887 : bessel_get_lim(double B, double L)
     171      320887 : { return maxss(2, L * exp(dbllambertW0(B))); }
     172             : 
     173          42 : static GEN vjbesselh(void* E, GEN z, long prec){return jbesselh((GEN)E,z,prec);}
     174         126 : static GEN vjbessel(void* E, GEN z, long prec) {return jbessel((GEN)E,z,prec);}
     175          42 : static GEN vibessel(void* E, GEN z, long prec) {return ibessel((GEN)E,z,prec);}
     176         126 : static GEN vnbessel(void* E, GEN z, long prec) {return nbessel((GEN)E,z,prec);}
     177          42 : static GEN vkbessel(void* E, GEN z, long prec) {return kbessel((GEN)E,z,prec);}
     178             : 
     179             : /* if J != 0 BesselJ, else BesselI. */
     180             : static GEN
     181      401245 : jbesselintern(GEN n, GEN z, long J, long prec)
     182             : {
     183      401245 :   const char *f = J? "besselj": "besseli";
     184             :   long i, ki;
     185      401245 :   pari_sp av = avma;
     186             :   GEN y;
     187             : 
     188      401245 :   switch(typ(z))
     189             :   {
     190      400839 :     case t_INT: case t_FRAC: case t_REAL: case t_COMPLEX:
     191             :     {
     192      400839 :       int flz0 = gequal0(z);
     193             :       long lim, k, precnew, bit;
     194             :       GEN p1, p2;
     195             :       double az, L;
     196             : 
     197      400839 :       i = precision(z); if (i) prec = i;
     198      400839 :       if (flz0 && gequal0(n)) return real_1(prec);
     199      400839 :       bit = prec2nbits(prec);
     200      400839 :       if (bessel_asymp(z, bit))
     201             :       {
     202       82758 :         GEN R = J? bessjasymp(n, z, bit): bessiasymp(n, z, bit);
     203       82758 :         if (typ(R) == t_COMPLEX && isexactzero(imag_i(n))
     204       82072 :                                 && gsigne(real_i(z)) > 0
     205       81918 :                                 && isexactzero(imag_i(z))) R = gcopy(gel(R,1));
     206       82758 :         return gerepileupto(av, R);
     207             :       }
     208      318081 :       p2 = gpow(gmul2n(z,-1),n,prec);
     209      318053 :       p2 = gdiv(p2, ggamma(gaddgs(n,1),prec));
     210      318053 :       if (flz0) return gerepileupto(av, p2);
     211      318053 :       az = dblmodulus(z); L = HALF_E * az;
     212      318053 :       precnew = prec;
     213      318053 :       if (az >= 1.0) precnew += 1 + nbits2extraprec((long)(az/M_LN2));
     214      318053 :       if (issmall(n,&ki)) {
     215      317136 :         k = labs(ki);
     216      317136 :         n = utoi(k);
     217             :       } else {
     218         833 :         i = precision(n);
     219         833 :         if (i && i < precnew) n = gtofp(n,precnew);
     220             :       }
     221      317969 :       z = gtofp(z,precnew);
     222      317969 :       lim = bessel_get_lim(prec2nbits_mul(prec,M_LN2/2) / L, L);
     223      317969 :       z = gmul2n(gsqr(z),-2); if (J) z = gneg(z);
     224      317969 :       p1 = gprec_wtrunc(_jbessel(n,z,lim), prec);
     225      317969 :       return gerepileupto(av, gmul(p2,p1));
     226             :     }
     227             : 
     228          14 :     case t_PADIC: pari_err_IMPL(stack_strcat("p-adic ",f));
     229         392 :     default:
     230             :     {
     231             :       long v, k, m;
     232         392 :       if (!(y = toser_i(z))) break;
     233         238 :       if (issmall(n,&ki)) n = utoi(labs(ki));
     234         210 :       y = gmul2n(gsqr(y),-2); if (J) y = gneg(y);
     235         210 :       v = valp(y);
     236         210 :       if (v < 0) pari_err_DOMAIN(f, "valuation", "<", gen_0, z);
     237         203 :       if (v == 0) pari_err_IMPL(stack_strcat(f, " around a!=0"));
     238         203 :       m = lg(y) - 2;
     239         203 :       k = m - (v >> 1);
     240         203 :       if (k <= 0) { set_avma(av); return scalarser(gen_1, varn(z), v); }
     241         203 :       setlg(y, k+2); return gerepileupto(av, _jbessel(n, y, m));
     242             :     }
     243             :   }
     244         154 :   return trans_evalgen(f, (void*)n, J? vjbessel: vibessel, z, prec);
     245             : }
     246             : GEN
     247      396885 : jbessel(GEN n, GEN z, long prec) { return jbesselintern(n,z,1,prec); }
     248             : GEN
     249         889 : ibessel(GEN n, GEN z, long prec) { return jbesselintern(n,z,0,prec); }
     250             : 
     251             : /* k > 0 */
     252             : static GEN
     253         119 : _jbesselh(long k, GEN z, long prec)
     254             : {
     255         119 :   GEN s, c, p0, p1, zinv = ginv(z);
     256             :   long i;
     257             : 
     258         119 :   gsincos(z,&s,&c,prec);
     259         119 :   p1 = gmul(zinv,s);
     260         119 :   p0 = p1; p1 = gmul(zinv,gsub(p0,c));
     261        1134 :   for (i = 2; i <= k; i++)
     262             :   {
     263        1015 :     GEN p2 = gsub(gmul(gmulsg(2*i-1,zinv), p1), p0);
     264        1015 :     p0 = p1; p1 = p2;
     265             :   }
     266         119 :   return p1;
     267             : }
     268             : 
     269             : /* J_{n+1/2}(z) */
     270             : GEN
     271         315 : jbesselh(GEN n, GEN z, long prec)
     272             : {
     273             :   long k, i;
     274             :   pari_sp av;
     275             :   GEN y;
     276             : 
     277         315 :   if (typ(n)!=t_INT) pari_err_TYPE("jbesselh",n);
     278         203 :   k = itos(n);
     279         203 :   if (k < 0) return jbessel(gadd(ghalf,n), z, prec);
     280             : 
     281         203 :   switch(typ(z))
     282             :   {
     283         133 :     case t_INT: case t_FRAC: case t_REAL: case t_COMPLEX:
     284             :     {
     285             :       long pr;
     286             :       GEN p1;
     287         133 :       if (gequal0(z))
     288             :       {
     289           7 :         av = avma;
     290           7 :         p1 = gmul(gsqrt(gdiv(z,mppi(prec)),prec),gpowgs(z,k));
     291           7 :         p1 = gdiv(p1, mulu_interval(k+1, 2*k+1)); /* x k! / (2k+1)! */
     292           7 :         return gerepileupto(av, gmul2n(p1,2*k));
     293             :       }
     294         126 :       if ( (pr = precision(z)) ) prec = pr;
     295         126 :       if (bessel_asymp(z, prec2nbits(prec)))
     296           7 :         return jbessel(gadd(ghalf,n), z, prec);
     297         119 :       y = cgetc(prec); av = avma;
     298         119 :       p1 = gsqrt(gdiv(z, Pi2n(-1,prec)), prec);
     299         119 :       if (!k)
     300          21 :         p1 = gmul(p1, gsinc(z, prec));
     301             :       else
     302             :       {
     303          98 :         long bits = BITS_IN_LONG + 2*k * (log2(k) -  dbllog2(z));
     304          98 :         if (bits > 0)
     305             :         {
     306          98 :           prec += nbits2extraprec(bits);
     307          98 :           if (pr) z = gtofp(z, prec);
     308             :         }
     309          98 :         p1 = gmul(p1, _jbesselh(k,z,prec));
     310             :       }
     311         119 :       set_avma(av); return affc_fixlg(p1, y);
     312             :     }
     313           0 :     case t_PADIC: pari_err_IMPL("p-adic jbesselh function");
     314          70 :     default:
     315             :     {
     316             :       long t, v;
     317          70 :       av = avma; if (!(y = toser_i(z))) break;
     318          35 :       if (gequal0(y)) return gerepileupto(av, gpowgs(y,k));
     319          35 :       v = valp(y);
     320          35 :       if (v < 0) pari_err_DOMAIN("besseljh","valuation", "<", gen_0, z);
     321          28 :       t = lg(y)-2;
     322          28 :       if (v) y = sertoser(y, t + (2*k+1)*v);
     323          28 :       if (!k)
     324           7 :         y = gsinc(y,prec);
     325             :       else
     326             :       {
     327          21 :         GEN T, a = _jbesselh(k, y, prec);
     328          21 :         if (v) y = sertoser(y, t + k*v); /* lower precision */
     329          21 :         y = gdiv(a, gpowgs(y, k));
     330          21 :         T = cgetg(k+1, t_VECSMALL);
     331         168 :         for (i = 1; i <= k; i++) T[i] = 2*i+1;
     332          21 :         y = gmul(y, zv_prod_Z(T));
     333             :       }
     334          28 :       return gerepileupto(av, y);
     335             :     }
     336             :   }
     337          35 :   return trans_evalgen("besseljh",(void*)n, vjbesselh, z, prec);
     338             : }
     339             : 
     340             : static GEN
     341           0 : kbessel2(GEN nu, GEN x, long prec)
     342             : {
     343           0 :   pari_sp av = avma;
     344           0 :   GEN p1, a, x2 = gshift(x,1);
     345             : 
     346           0 :   a = gtofp(gaddgs(gshift(nu,1), 1), prec);
     347           0 :   p1 = hyperu(gshift(a,-1), a, x2, prec);
     348           0 :   p1 = gmul(gmul(p1, gpow(x2,nu,prec)), sqrtr(mppi(prec)));
     349           0 :   return gerepileupto(av, gmul(p1, gexp(gneg(x),prec)));
     350             : }
     351             : 
     352             : /* special case of hyperu */
     353             : static GEN
     354          14 : kbessel1(GEN nu, GEN gx, long prec)
     355             : {
     356             :   GEN x, y, zf, r, u, pi, nu2;
     357          14 :   long bit, k, k2, n2, n, l = (typ(gx)==t_REAL)? realprec(gx): prec;
     358             :   pari_sp av;
     359             : 
     360          14 :   if (typ(nu)==t_COMPLEX) return kbessel2(nu, gx, l);
     361          14 :   y = cgetr(l); av = avma;
     362          14 :   x = gtofp(gx, l);
     363          14 :   nu = gtofp(nu,l); nu2 = sqrr(nu);
     364          14 :   shiftr_inplace(nu2,2); togglesign(nu2); /* nu2 = -4nu^2 */
     365          14 :   n = (long) (prec2nbits_mul(l,M_LN2) + M_PI*fabs(rtodbl(nu))) / 2;
     366          14 :   bit = prec2nbits(l) - 1;
     367          14 :   l += EXTRAPRECWORD;
     368          14 :   pi = mppi(l); n2 = n<<1; r = gmul2n(x,1);
     369          14 :   if (cmprs(x, n) < 0)
     370             :   {
     371          14 :     pari_sp av2 = avma;
     372          14 :     GEN q, v, c, s = real_1(l), t = real_0(l);
     373        1246 :     for (k = n2, k2 = 2*n2-1; k > 0; k--, k2 -= 2)
     374             :     {
     375        1232 :       GEN ak = divri(addri(nu2, sqru(k2)), mulss(n2<<2, -k));
     376        1232 :       s = addsr(1, mulrr(ak,s));
     377        1232 :       t = addsr(k2,mulrr(ak,t));
     378        1232 :       if (gc_needed(av2,3)) gerepileall(av2, 2, &s,&t);
     379             :     }
     380          14 :     shiftr_inplace(t, -1);
     381          14 :     q = utor(n2, l);
     382          14 :     zf = sqrtr(divru(pi,n2));
     383          14 :     u = gprec_wensure(mulrr(zf, s), l);
     384          14 :     v = gprec_wensure(divrs(addrr(mulrr(t,zf),mulrr(u,nu)),-n2), l);
     385             :     for(;;)
     386         301 :     {
     387         315 :       GEN p1, e, f, d = real_1(l);
     388             :       pari_sp av3;
     389         315 :       c = divur(5,q); if (expo(c) >= -1) c = real2n(-1,l);
     390         315 :       p1 = subsr(1, divrr(r,q)); if (cmprr(c,p1)>0) c = p1;
     391         315 :       togglesign(c); av3 = avma;
     392         315 :       e = u;
     393         315 :       f = v;
     394         315 :       for (k = 1;; k++)
     395       34475 :       {
     396       34790 :         GEN w = addrr(gmul2n(mulur(2*k-1,u), -1), mulrr(subrs(q,k),v));
     397       34790 :         w = addrr(w, mulrr(nu, subrr(u,gmul2n(v,1))));
     398       34790 :         u = divru(mulrr(q,v), k);
     399       34790 :         v = divru(w,k);
     400       34790 :         d = mulrr(d,c);
     401       34790 :         e = addrr(e, mulrr(d,u));
     402       34790 :         f = addrr(f, p1 = mulrr(d,v));
     403       34790 :         if (expo(p1) - expo(f) <= 1-prec2nbits(realprec(p1))) break;
     404       34475 :         if (gc_needed(av3,3)) gerepileall(av3,5,&u,&v,&d,&e,&f);
     405             :       }
     406         315 :       u = e;
     407         315 :       v = f;
     408         315 :       q = mulrr(q, addrs(c,1));
     409         315 :       if (expo(r) - expo(subrr(q,r)) >= bit) break;
     410         301 :       gerepileall(av2, 3, &u,&v,&q);
     411             :     }
     412          14 :     u = mulrr(u, gpow(divru(x,n),nu,prec));
     413             :   }
     414             :   else
     415             :   {
     416           0 :     GEN s, zz = ginv(gmul2n(r,2));
     417           0 :     pari_sp av2 = avma;
     418           0 :     s = real_1(l);
     419           0 :     for (k = n2, k2 = 2*n2-1; k > 0; k--, k2 -= 2)
     420             :     {
     421           0 :       GEN ak = divru(mulrr(addri(nu2, sqru(k2)), zz), k);
     422           0 :       s = subsr(1, mulrr(ak,s));
     423           0 :       if (gc_needed(av2,3)) s = gerepileuptoleaf(av2, s);
     424             :     }
     425           0 :     zf = sqrtr(divrr(pi,r));
     426           0 :     u = mulrr(s, zf);
     427             :   }
     428          14 :   affrr(mulrr(u, mpexp(negr(x))), y);
     429          14 :   set_avma(av); return y;
     430             : }
     431             : 
     432             : /*   sum_{k=0}^m x^k (H(k)+H(k+n)) / (k! (k+n)!)
     433             :  * + sum_{k=0}^{n-1} (-x)^(k-n) (n-k-1)!/k! */
     434             : static GEN
     435        3002 : _kbessel(long n, GEN x, long m, long prec)
     436             : {
     437             :   GEN p1, p2, s, H;
     438        3002 :   long k, M = m + n, exact = (M <= bit_accuracy(prec));
     439             :   pari_sp av;
     440             : 
     441        3002 :   H = cgetg(M+2,t_VEC); gel(H,1) = gen_0;
     442        3002 :   if (exact)
     443             :   {
     444        3002 :     gel(H,2) = s = gen_1;
     445      107856 :     for (k=2; k<=M; k++) gel(H,k+1) = s = gdivgu(gaddsg(1,gmulsg(k,s)),k);
     446             :   }
     447             :   else
     448             :   {
     449           0 :     gel(H,2) = s = real_1(prec);
     450           0 :     for (k=2; k<=M; k++) gel(H,k+1) = s = divru(addsr(1,mulur(k,s)),k);
     451             :   }
     452        3002 :   s = gadd(gel(H,m+1), gel(H,M+1)); av = avma;
     453      109378 :   for (k = m; k > 0; k--)
     454             :   {
     455      106376 :     s = gadd(gadd(gel(H,k),gel(H,k+n)), gdiv(gmul(x,s),mulss(k,k+n)));
     456      106376 :     if (gc_needed(av,1))
     457             :     {
     458           0 :       if (DEBUGMEM>1) pari_warn(warnmem,"_kbessel");
     459           0 :       s = gerepileupto(av, s);
     460             :     }
     461             :   }
     462        3002 :   p1 = exact? mpfact(n): mpfactr(n,prec);
     463        3002 :   s = gdiv(s,p1);
     464        3002 :   if (n)
     465             :   {
     466         472 :     x = gneg(ginv(x));
     467         472 :     p2 = gmulsg(n, gdiv(x,p1));
     468         472 :     s = gadd(s,p2);
     469        1480 :     for (k=n-1; k>0; k--)
     470             :     {
     471        1008 :       p2 = gmul(p2, gmul(mulss(k,n-k),x));
     472        1008 :       s = gadd(s,p2);
     473             :     }
     474             :   }
     475        3002 :   return s;
     476             : }
     477             : 
     478             : /* N = 1: Bessel N, else Bessel K */
     479             : static GEN
     480        4578 : kbesselintern(GEN n, GEN z, long N, long prec)
     481             : {
     482        4578 :   const char *f = N? "besseln": "besselk";
     483             :   long i, k, ki, lim, precnew, fl2, ex, bit;
     484        4578 :   pari_sp av = avma;
     485             :   GEN p1, p2, y, p3, pp, pm, s, c;
     486             :   double az;
     487             : 
     488        4578 :   switch(typ(z))
     489             :   {
     490        4193 :     case t_INT: case t_FRAC: case t_REAL: case t_COMPLEX:
     491        4193 :       if (gequal0(z)) pari_err_DOMAIN(f, "argument", "=", gen_0, z);
     492        4193 :       i = precision(z); if (i) prec = i;
     493        4193 :       i = precision(n); if (i && prec > i) prec = i;
     494        4193 :       bit = prec2nbits(prec);
     495        4193 :       if (bessel_asymp(z, bit))
     496             :       {
     497         932 :         GEN R = N? bessyasymp(n, z, bit): besskasymp(n, z, bit);
     498         932 :         if (typ(R) == t_COMPLEX && isexactzero(imag_i(n))
     499         225 :                                 && gsigne(real_i(z)) > 0
     500          85 :                                 && isexactzero(imag_i(z))) R = gcopy(gel(R,1));
     501         932 :         return gerepileupto(av, R);
     502             :       }
     503             :       /* heuristic threshold */
     504        3261 :       if (!N && !gequal0(n) && gexpo(z) > bit/16 + gexpo(n))
     505          14 :         return kbessel1(n,z,prec);
     506        3240 :       az = dblmodulus(z); precnew = prec;
     507        3240 :       if (az >= 1) precnew += 1 + nbits2extraprec((long)((N?az:2*az)/M_LN2));
     508        3240 :       z = gtofp(z, precnew);
     509        3240 :       if (issmall(n,&ki))
     510             :       {
     511        2918 :         GEN z2 = gmul2n(z, -1), Z;
     512        2918 :         double B, L = HALF_E * az;
     513        2918 :         k = labs(ki);
     514        2918 :         B = prec2nbits_mul(prec,M_LN2/2) / L;
     515        2918 :         if (!N) B += 0.367879; /* exp(-1) */
     516        2918 :         lim = bessel_get_lim(B, L);
     517        2918 :         Z = gsqr(z2); if (N) Z = gneg(Z);
     518        2918 :         p1 = gmul(gpowgs(z2,k), _kbessel(k, Z, lim, precnew));
     519        2918 :         p2 = gadd(mpeuler(precnew), glog(z2,precnew));
     520        2918 :         p3 = jbesselintern(stoi(k),z,N,precnew);
     521        2918 :         p2 = gsub(gmul2n(p1,-1),gmul(p2,p3));
     522        2918 :         p2 = gprec_wtrunc(p2, prec);
     523        2918 :         if (N)
     524             :         {
     525         510 :           p2 = gdiv(p2, Pi2n(-1,prec));
     526         510 :           if (ki >= 0 || !odd(k)) p2 = gneg(p2);
     527             :         } else
     528        2408 :           if (odd(k)) p2 = gneg(p2);
     529        2918 :         return gerepilecopy(av, p2);
     530             :       }
     531             : 
     532         259 :       n = gtofp(n, precnew);
     533         259 :       gsincos(gmul(n,mppi(precnew)), &s,&c,precnew);
     534         259 :       ex = gexpo(s);
     535         259 :       if (ex < 0) precnew += nbits2extraprec(N? -ex: -2*ex);
     536         259 :       if (i && i < precnew) {
     537          84 :         n = gtofp(n,precnew);
     538          84 :         z = gtofp(z,precnew);
     539          84 :         gsincos(gmul(n,mppi(precnew)), &s,&c,precnew);
     540             :       }
     541             : 
     542         259 :       pp = jbesselintern(n,      z,N,precnew);
     543         259 :       pm = jbesselintern(gneg(n),z,N,precnew);
     544         259 :       if (N)
     545         189 :         p1 = gsub(gmul(c,pp),pm);
     546             :       else
     547          70 :         p1 = gmul(gsub(pm,pp), Pi2n(-1,precnew));
     548         259 :       p1 = gdiv(p1, s);
     549         259 :       return gerepilecopy(av, gprec_wtrunc(p1,prec));
     550             : 
     551          14 :     case t_PADIC: pari_err_IMPL(stack_strcat("p-adic ",f));
     552         371 :     default:
     553         371 :       if (!(y = toser_i(z))) break;
     554         217 :       if (issmall(n,&ki))
     555             :       {
     556         105 :         long v, mv, k = labs(ki), m = lg(y)-2;
     557         105 :         y = gmul2n(gsqr(y),-2); if (N) y = gneg(y);
     558         105 :         v = valp(y);
     559         105 :         if (v < 0) pari_err_DOMAIN(f, "valuation", "<", gen_0, z);
     560          91 :         if (v == 0) pari_err_IMPL(stack_strcat(f, " around a!=0"));
     561          91 :         mv = m - (v >> 1);
     562          91 :         if (mv <= 0) { set_avma(av); return scalarser(gen_1, varn(z), v); }
     563          84 :         setlg(y, mv+2); return gerepilecopy(av, _kbessel(k, y, m, prec));
     564             :       }
     565          98 :       if (!issmall(gmul2n(n,1),&ki))
     566          70 :         pari_err_DOMAIN(f, "2n mod Z", "!=", gen_0, n);
     567          28 :       k = labs(ki); n = gmul2n(stoi(k),-1);
     568          28 :       fl2 = (k&3)==1;
     569          28 :       pm = jbesselintern(gneg(n), y, N, prec);
     570          28 :       if (N) p1 = pm;
     571             :       else
     572             :       {
     573           7 :         pp = jbesselintern(n, y, N, prec);
     574           7 :         p2 = gpowgs(y,-k); if (fl2 == 0) p2 = gneg(p2);
     575           7 :         p3 = gmul2n(diviiexact(mpfact(k + 1),mpfact((k + 1) >> 1)),-(k + 1));
     576           7 :         p3 = gdivgu(gmul2n(gsqr(p3),1),k);
     577           7 :         p2 = gmul(p2,p3);
     578           7 :         p1 = gsub(pp,gmul(p2,pm));
     579             :       }
     580          28 :       return gerepileupto(av, fl2? gneg(p1): gcopy(p1));
     581             :   }
     582         154 :   return trans_evalgen(f, (void*)n, N? vnbessel: vkbessel, z, prec);
     583             : }
     584             : 
     585             : GEN
     586        3108 : kbessel(GEN n, GEN z, long prec) { return kbesselintern(n,z,0,prec); }
     587             : GEN
     588        1470 : ybessel(GEN n, GEN z, long prec) { return kbesselintern(n,z,1,prec); }
     589             : /* J + iN */
     590             : GEN
     591         224 : hbessel1(GEN n, GEN z, long prec)
     592             : {
     593         224 :   pari_sp av = avma;
     594         224 :   GEN J = jbessel(n,z,prec);
     595         196 :   GEN Y = ybessel(n,z,prec);
     596         182 :   return gerepileupto(av, gadd(J, mulcxI(Y)));
     597             : }
     598             : /* J - iN */
     599             : GEN
     600         224 : hbessel2(GEN n, GEN z, long prec)
     601             : {
     602         224 :   pari_sp av = avma;
     603         224 :   GEN J = jbessel(n,z,prec);
     604         196 :   GEN Y = ybessel(n,z,prec);
     605         182 :   return gerepileupto(av, gadd(J, mulcxmI(Y)));
     606             : }
     607             : 
     608             : static GEN
     609          28 : besselrefine(GEN z, GEN nu, GEN (*B)(GEN,GEN,long), long bit)
     610             : {
     611          28 :   GEN z0 = gprec_w(z, DEFAULTPREC), nu1 = gaddgs(nu, 1), t;
     612          28 :   long e, n, c, j, prec = DEFAULTPREC;
     613             : 
     614          28 :   t = gdiv(B(nu1, z0, prec), B(nu, z0, prec));
     615          28 :   t = gadd(z0, gdiv(gsub(gsqr(z0), gsqr(nu)), gsub(gdiv(nu, z0), t)));
     616          28 :   e = gexpo(t) - 2 * gexpo(z0) - 1; if (e < 0) e = 0;
     617          28 :   n = expu(bit + 32 - e);
     618          28 :   c = 1 + e + ((bit - e) >> n);
     619         224 :   for (j = 1; j <= n; j++)
     620             :   {
     621         196 :     c = 2 * c - e;
     622         196 :     prec = nbits2prec(c); z = gprec_w(z, prec);
     623         196 :     t = gdiv(B(nu1, z, prec), B(nu, z, prec));
     624         196 :     z = gsub(z, ginv(gsub(gdiv(nu, z), t)));
     625             :   }
     626          28 :   return gprec_w(z, nbits2prec(bit));
     627             : }
     628             : 
     629             : static GEN
     630          42 : besselzero(GEN nu, long n, GEN (*B)(GEN,GEN,long), long bit)
     631             : {
     632          42 :   pari_sp av = avma;
     633          42 :   long prec = nbits2prec(bit), a;
     634             :   GEN z;
     635          42 :   if (n <= 0) pari_err_DOMAIN("besselzero", "n", "<=", gen_0, stoi(n));
     636          28 :   if (n > LONG_MAX / 4) pari_err_OVERFLOW("besselzero");
     637          28 :   a = 4 * n - (B == jbessel? 1: 3);
     638          28 :   z = gmul(mppi(prec), gmul2n(gaddgs(gmul2n(nu,1), a), -2));
     639          28 :   return gerepilecopy(av, besselrefine(z, nu, B, bit));
     640             : }
     641             : GEN
     642          21 : besseljzero(GEN nu, long k, long b) { return besselzero(nu, k, jbessel, b); }
     643             : GEN
     644          21 : besselyzero(GEN nu, long k, long b) { return besselzero(nu, k, ybessel, b); }
     645             : 
     646             : /***********************************************************************/
     647             : /**                    INCOMPLETE GAMMA FUNCTION                      **/
     648             : /***********************************************************************/
     649             : /* mx ~ |x|, b = bit accuracy */
     650             : static int
     651       16422 : gamma_use_asymp(GEN x, long b)
     652             : {
     653             :   long e;
     654       16422 :   if (is_real_t(typ(x)))
     655             :   {
     656       12390 :     pari_sp av = avma;
     657       12390 :     return gc_int(av, gcmpgs(R_abs_shallow(x), 3*b / 4) >= 0);
     658             :   }
     659        4032 :   e = gexpo(x); return e >= b || dblmodulus(x) >= 3*b / 4;
     660             : }
     661             : /* x a t_REAL */
     662             : static GEN
     663          28 : eint1r_asymp(GEN x, GEN expx, long prec)
     664             : {
     665          28 :   pari_sp av = avma, av2;
     666             :   GEN S, q, z, ix;
     667          28 :   long oldeq = LONG_MAX, esx = -prec2nbits(prec), j;
     668             : 
     669          28 :   if (realprec(x) < prec + EXTRAPREC64) x = rtor(x, prec+EXTRAPREC64);
     670          28 :   ix = invr(x); q = z = negr(ix);
     671          28 :   av2 = avma; S = addrs(q, 1);
     672          28 :   for (j = 2;; j++)
     673        1211 :   {
     674        1239 :     long eq = expo(q); if (eq < esx) break;
     675        1211 :     if ((j & 3) == 0)
     676             :     { /* guard against divergence */
     677         294 :       if (eq > oldeq) return gc_NULL(av); /* regressing, abort */
     678         294 :       oldeq = eq;
     679             :     }
     680        1211 :     q = mulrr(q, mulru(z, j)); S = addrr(S, q);
     681        1211 :     if (gc_needed(av2, 1)) gerepileall(av2, 2, &S, &q);
     682             :   }
     683          28 :   if (DEBUGLEVEL > 2) err_printf("eint1: using asymp\n");
     684          28 :   S = expx? divrr(S, expx): mulrr(S, mpexp(negr(x)));
     685          28 :   return gerepileuptoleaf(av, mulrr(S, ix));
     686             : }
     687             : /* cf incgam_asymp(0, x); z = -1/x
     688             :  *   exp(-x)/x * (1 + z + 2! z^2 + ...) */
     689             : static GEN
     690         105 : eint1_asymp(GEN x, GEN expx, long prec)
     691             : {
     692         105 :   pari_sp av = avma, av2;
     693             :   GEN S, q, z, ix;
     694         105 :   long oldeq = LONG_MAX, esx = -prec2nbits(prec), j;
     695             : 
     696         105 :   if (typ(x) != t_REAL) x = gtofp(x, prec+EXTRAPREC64);
     697         105 :   if (typ(x) == t_REAL) return eint1r_asymp(x, expx, prec);
     698         105 :   ix = ginv(x); q = z = gneg_i(ix);
     699         105 :   av2 = avma; S = gaddgs(q, 1);
     700         105 :   for (j = 2;; j++)
     701        5824 :   {
     702        5929 :     long eq = gexpo(q); if (eq < esx) break;
     703        5824 :     if ((j & 3) == 0)
     704             :     { /* guard against divergence */
     705        1442 :       if (eq > oldeq) return gc_NULL(av); /* regressing, abort */
     706        1442 :       oldeq = eq;
     707             :     }
     708        5824 :     q = gmul(q, gmulgu(z, j)); S = gadd(S, q);
     709        5824 :     if (gc_needed(av2, 1)) gerepileall(av2, 2, &S, &q);
     710             :   }
     711         105 :   if (DEBUGLEVEL > 2) err_printf("eint1: using asymp\n");
     712         105 :   S = expx? gdiv(S, expx): gmul(S, gexp(gneg_i(x), prec));
     713         105 :   return gerepileupto(av, gmul(S, ix));
     714             : }
     715             : 
     716             : /* eint1(x) = incgam(0, x); typ(x) = t_REAL, x > 0 */
     717             : static GEN
     718        6188 : eint1p(GEN x, GEN expx)
     719             : {
     720             :   pari_sp av;
     721        6188 :   long l = realprec(x), bit = prec2nbits(l), prec, i;
     722             :   double mx;
     723             :   GEN z, S, t, H, run;
     724             : 
     725        6188 :   if (gamma_use_asymp(x, bit)
     726          28 :       && (z = eint1r_asymp(x, expx, l))) return z;
     727        6160 :   mx = rtodbl(x);
     728        6160 :   prec = l + nbits2extraprec((mx+log(mx))/M_LN2 + 10);
     729        6160 :   bit = prec2nbits(prec);
     730        6160 :   run = real_1(prec); x = rtor(x, prec);
     731        6160 :   av = avma; S = z = t = H = run;
     732      623624 :   for (i = 2; expo(S) - expo(t) <= bit; i++)
     733             :   {
     734      617464 :     H = addrr(H, divru(run,i)); /* H = sum_{k<=i} 1/k */
     735      617464 :     z = divru(mulrr(x,z), i);   /* z = x^(i-1)/i! */
     736      617464 :     t = mulrr(z, H); S = addrr(S, t);
     737      617464 :     if ((i & 0x1ff) == 0) gerepileall(av, 4, &z,&t,&S,&H);
     738             :   }
     739        6160 :   return subrr(mulrr(x, divrr(S,expx? expx: mpexp(x))),
     740             :                addrr(mplog(x), mpeuler(prec)));
     741             : }
     742             : /* eint1(x) = incgam(0, x); typ(x) = t_REAL, x < 0
     743             :  * rewritten from code contributed by Manfred Radimersky */
     744             : static GEN
     745         140 : eint1m(GEN x, GEN expx)
     746             : {
     747         140 :   GEN p1, q, S, y, z = cgetg(3, t_COMPLEX);
     748         140 :   long l  = realprec(x), n  = prec2nbits(l), j;
     749         140 :   pari_sp av = avma;
     750             : 
     751         140 :   y  = rtor(x, l + EXTRAPREC64); setsigne(y,1); /* |x| */
     752         140 :   if (gamma_use_asymp(y, n))
     753             :   { /* ~eint1_asymp: asymptotic expansion */
     754          14 :     p1 = q = invr(y); S = addrs(q, 1);
     755         560 :     for (j = 2; expo(q) >= -n; j++) {
     756         546 :       q = mulrr(q, mulru(p1, j));
     757         546 :       S = addrr(S, q);
     758             :     }
     759          14 :     y  = mulrr(p1, expx? divrr(S, expx): mulrr(S, mpexp(y)));
     760             :   }
     761             :   else
     762             :   {
     763         126 :     p1 = q = S = y;
     764       24248 :     for (j = 2; expo(q) - expo(S) >= -n; j++) {
     765       24122 :       p1 = mulrr(y, divru(p1, j)); /* (-x)^j/j! */
     766       24122 :       q = divru(p1, j);
     767       24122 :       S = addrr(S, q);
     768             :     }
     769         126 :     y  = addrr(S, addrr(logr_abs(x), mpeuler(l)));
     770             :   }
     771         140 :   y = gerepileuptoleaf(av, y); togglesign(y);
     772         140 :   gel(z, 1) = y;
     773         140 :   y = mppi(l); setsigne(y, -1);
     774         140 :   gel(z, 2) = y; return z;
     775             : }
     776             : 
     777             : /* real(z*log(z)-z), z = x+iy */
     778             : static double
     779        7714 : mygamma(double x, double y)
     780             : {
     781        7714 :   if (x == 0.) return -(M_PI/2)*fabs(y);
     782        7714 :   return (x/2)*log(x*x+y*y)-x-y*atan(y/x);
     783             : }
     784             : 
     785             : /* x^s exp(-x) */
     786             : static GEN
     787       10066 : expmx_xs(GEN s, GEN x, GEN logx, long prec)
     788             : {
     789             :   GEN z;
     790       10066 :   long ts = typ(s);
     791       10066 :   if (ts == t_INT || (ts == t_FRAC && absequaliu(gel(s,2), 2)))
     792        5257 :     z = gmul(gexp(gneg(x), prec), gpow(x, s, prec));
     793             :   else
     794        4809 :     z = gexp(gsub(gmul(s, logx? logx: glog(x,prec+EXTRAPREC64)), x), prec);
     795       10066 :   return z;
     796             : }
     797             : 
     798             : /* Not yet: doesn't work at low accuracy
     799             : #define INCGAM_CF
     800             : */
     801             : 
     802             : #ifdef INCGAM_CF
     803             : /* Is s very close to a nonpositive integer ? */
     804             : static int
     805             : isgammapole(GEN s, long bitprec)
     806             : {
     807             :   pari_sp av = avma;
     808             :   GEN t = imag_i(s);
     809             :   long e, b = bitprec - 10;
     810             : 
     811             :   if (gexpo(t) > - b) return 0;
     812             :   s = real_i(s);
     813             :   if (gsigne(s) > 0 && gexpo(s) > -b) return 0;
     814             :   (void)grndtoi(s, &e); return gc_bool(av, e < -b);
     815             : }
     816             : 
     817             : /* incgam using the continued fraction. x a t_REAL or t_COMPLEX, mx ~ |x|.
     818             :  * Assume precision(s), precision(x) >= prec */
     819             : static GEN
     820             : incgam_cf(GEN s, GEN x, double mx, long prec)
     821             : {
     822             :   GEN ms, y, S;
     823             :   long n, i, j, LS, bitprec = prec2nbits(prec);
     824             :   double rs, is, m;
     825             : 
     826             :   if (typ(s) == t_COMPLEX)
     827             :   {
     828             :     rs = gtodouble(gel(s,1));
     829             :     is = gtodouble(gel(s,2));
     830             :   }
     831             :   else
     832             :   {
     833             :     rs = gtodouble(s);
     834             :     is = 0.;
     835             :   }
     836             :   if (isgammapole(s, bitprec)) LS = 0;
     837             :   else
     838             :   {
     839             :     double bit,  LGS = mygamma(rs,is);
     840             :     LS = LGS <= 0 ? 0: ceil(LGS);
     841             :     bit = (LGS - (rs-1)*log(mx) + mx)/M_LN2;
     842             :     if (bit > 0)
     843             :     {
     844             :       prec += nbits2extraprec((long)bit);
     845             :       x = gtofp(x, prec);
     846             :       if (isinexactreal(s)) s = gtofp(s, prec);
     847             :     }
     848             :   }
     849             :   /* |ln(2*gamma(s)*sin(s*Pi))| <= ln(2) + |lngamma(s)| + |Im(s)*Pi|*/
     850             :   m = bitprec*M_LN2 + LS + M_LN2 + fabs(is)*M_PI + mx;
     851             :   if (rs < 1) m += (1 - rs)*log(mx);
     852             :   m /= 4;
     853             :   n = (long)(1 + m*m/mx);
     854             :   y = expmx_xs(gsubgs(s,1), x, NULL, prec);
     855             :   if (rs >= 0 && bitprec >= 512)
     856             :   {
     857             :     GEN A = cgetg(n+1, t_VEC), B = cgetg(n+1, t_VEC);
     858             :     ms = gsubsg(1, s);
     859             :     for (j = 1; j <= n; ++j)
     860             :     {
     861             :       gel(A,j) = ms;
     862             :       gel(B,j) = gmulsg(j, gsubgs(s,j));
     863             :       ms = gaddgs(ms, 2);
     864             :     }
     865             :     S = contfraceval_inv(mkvec2(A,B), x, -1);
     866             :   }
     867             :   else
     868             :   {
     869             :     GEN x_s = gsub(x, s);
     870             :     pari_sp av2 = avma;
     871             :     S = gdiv(gsubgs(s,n), gaddgs(x_s,n<<1));
     872             :     for (i=n-1; i >= 1; i--)
     873             :     {
     874             :       S = gdiv(gsubgs(s,i), gadd(gaddgs(x_s,i<<1),gmulsg(i,S)));
     875             :       if (gc_needed(av2,3))
     876             :       {
     877             :         if(DEBUGMEM>1) pari_warn(warnmem,"incgam_cf");
     878             :         S = gerepileupto(av2, S);
     879             :       }
     880             :     }
     881             :     S = gaddgs(S,1);
     882             :   }
     883             :   return gmul(y, S);
     884             : }
     885             : #endif
     886             : 
     887             : static double
     888        5642 : findextraincgam(GEN s, GEN x)
     889             : {
     890        5642 :   double sig = gtodouble(real_i(s)), t = gtodouble(imag_i(s));
     891        5642 :   double xr = gtodouble(real_i(x)), xi = gtodouble(imag_i(x));
     892        5642 :   double exd = 0., Nx = xr*xr + xi*xi, D = Nx - t*t;
     893             :   long n;
     894             : 
     895        5642 :   if (xr < 0)
     896             :   {
     897         833 :     long ex = gexpo(x);
     898         833 :     if (ex > 0 && ex > gexpo(s)) exd = sqrt(Nx)*log(Nx)/2; /* |x| log |x| */
     899             :   }
     900        5642 :   if (D <= 0.) return exd;
     901        4200 :   n = (long)(sqrt(D)-sig);
     902        4200 :   if (n <= 0) return exd;
     903        1512 :   return maxdd(exd, (n*log(Nx)/2 - mygamma(sig+n, t) + mygamma(sig, t)) / M_LN2);
     904             : }
     905             : 
     906             : /* use exp(-x) * (x^s/s) * sum_{k >= 0} x^k / prod(i=1, k, s+i) */
     907             : static GEN
     908        5649 : incgamc_i(GEN s, GEN x, long *ptexd, long prec)
     909             : {
     910             :   GEN S, t, y;
     911             :   long l, n, i, exd;
     912        5649 :   pari_sp av = avma, av2;
     913             : 
     914        5649 :   if (gequal0(x))
     915             :   {
     916           7 :     if (ptexd) *ptexd = 0.;
     917           7 :     return gtofp(x, prec);
     918             :   }
     919        5642 :   l = precision(x);
     920        5642 :   if (!l) l = prec;
     921        5642 :   n = -prec2nbits(l)-1;
     922        5642 :   exd = (long)findextraincgam(s, x);
     923        5642 :   if (ptexd) *ptexd = exd;
     924        5642 :   if (exd > 0)
     925             :   {
     926        1393 :     long p = l + nbits2extraprec(exd);
     927        1393 :     x = gtofp(x, p);
     928        1393 :     if (isinexactreal(s)) s = gtofp(s, p);
     929             :   }
     930        4249 :   else x = gtofp(x, l+EXTRAPREC64);
     931        5642 :   av2 = avma;
     932        5642 :   S = gdiv(x, gaddsg(1,s));
     933        5642 :   t = gaddsg(1, S);
     934      753599 :   for (i=2; gexpo(S) >= n; i++)
     935             :   {
     936      747957 :     S = gdiv(gmul(x,S), gaddsg(i,s)); /* x^i / ((s+1)...(s+i)) */
     937      747957 :     t = gadd(S,t);
     938      747957 :     if (gc_needed(av2,3))
     939             :     {
     940           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"incgamc");
     941           0 :       gerepileall(av2, 2, &S, &t);
     942             :     }
     943             :   }
     944        5642 :   y = expmx_xs(s, x, NULL, prec);
     945        5642 :   return gerepileupto(av, gmul(gdiv(y,s), t));
     946             : }
     947             : 
     948             : GEN
     949        1456 : incgamc(GEN s, GEN x, long prec)
     950        1456 : { return incgamc_i(s, x, NULL, prec); }
     951             : 
     952             : /* incgamma using asymptotic expansion:
     953             :  *   exp(-x)x^(s-1)(1 + (s-1)/x + (s-1)(s-2)/x^2 + ...) */
     954             : static GEN
     955        2716 : incgam_asymp(GEN s, GEN x, long prec)
     956             : {
     957        2716 :   pari_sp av = avma, av2;
     958             :   GEN S, q, cox, invx;
     959        2716 :   long oldeq = LONG_MAX, eq, esx, j;
     960        2716 :   int flint = (typ(s) == t_INT && signe(s) > 0);
     961             : 
     962        2716 :   x = gtofp(x,prec+EXTRAPREC64);
     963        2716 :   invx = ginv(x);
     964        2716 :   esx = -prec2nbits(prec);
     965        2716 :   av2 = avma;
     966        2716 :   q = gmul(gsubgs(s,1), invx);
     967        2716 :   S = gaddgs(q, 1);
     968        2716 :   for (j = 2;; j++)
     969             :   {
     970      123879 :     eq = gexpo(q); if (eq < esx) break;
     971      121282 :     if (!flint && (j & 3) == 0)
     972             :     { /* guard against divergence */
     973       15778 :       if (eq > oldeq) return gc_NULL(av); /* regressing, abort */
     974       15659 :       oldeq = eq;
     975             :     }
     976      121163 :     q = gmul(q, gmul(gsubgs(s,j), invx));
     977      121163 :     S = gadd(S, q);
     978      121163 :     if (gc_needed(av2, 1)) gerepileall(av2, 2, &S, &q);
     979             :   }
     980        2597 :   if (DEBUGLEVEL > 2) err_printf("incgam: using asymp\n");
     981        2597 :   cox = expmx_xs(gsubgs(s,1), x, NULL, prec);
     982        2597 :   return gerepileupto(av, gmul(cox, S));
     983             : }
     984             : 
     985             : /* gasx = incgam(s-n,x). Compute incgam(s,x)
     986             :  * = (s-1)(s-2)...(s-n)gasx + exp(-x)x^(s-1) *
     987             :  *   (1 + (s-1)/x + ... + (s-1)(s-2)...(s-n+1)/x^(n-1)) */
     988             : static GEN
     989         546 : incgam_asymp_partial(GEN s, GEN x, GEN gasx, long n, long prec)
     990             : {
     991             :   pari_sp av;
     992         546 :   GEN S, q, cox, invx, s1 = gsubgs(s, 1), sprod;
     993             :   long j;
     994         546 :   cox = expmx_xs(s1, x, NULL, prec);
     995         546 :   if (n == 1) return gadd(cox, gmul(s1, gasx));
     996         546 :   invx = ginv(x);
     997         546 :   av = avma;
     998         546 :   q = gmul(s1, invx);
     999         546 :   S = gaddgs(q, 1);
    1000       52164 :   for (j = 2; j < n; j++)
    1001             :   {
    1002       51618 :     q = gmul(q, gmul(gsubgs(s, j), invx));
    1003       51618 :     S = gadd(S, q);
    1004       51618 :     if (gc_needed(av, 2)) gerepileall(av, 2, &S, &q);
    1005             :   }
    1006         546 :   sprod = gmul(gmul(q, gpowgs(x, n-1)), gsubgs(s, n));
    1007         546 :   return gadd(gmul(cox, S), gmul(sprod, gasx));
    1008             : }
    1009             : 
    1010             : /* Assume s != 0; called when Re(s) <= 1/2 */
    1011             : static GEN
    1012        2401 : incgamspec(GEN s, GEN x, GEN g, long prec)
    1013             : {
    1014        2401 :   GEN q, S, cox = gen_0, P, sk, S1, S2, S3, F3, logx, mx;
    1015        2401 :   long n, esk, E, k = itos(ground(gneg(real_i(s)))); /* >= 0 */
    1016             : 
    1017        2401 :   if (k && gexpo(x) > 0)
    1018             :   {
    1019         245 :     GEN xk = gdivgu(x, k);
    1020         245 :     long bitprec = prec2nbits(prec);
    1021         245 :     double d = (gexpo(xk) > bitprec)? bitprec*M_LN2: log(dblmodulus(xk));
    1022         245 :     d = k * (d + 1) / M_LN2;
    1023         245 :     if (d > 0) prec += nbits2extraprec((long)d);
    1024         245 :     if (isinexactreal(s)) s = gtofp(s, prec);
    1025             :   }
    1026        2401 :   x = gtofp(x, maxss(precision(x), prec) + EXTRAPREC64);
    1027        2401 :   sk = gaddgs(s, k); /* |Re(sk)| <= 1/2 */
    1028        2401 :   logx = glog(x, prec);
    1029        2401 :   mx = gneg(x);
    1030        2401 :   if (k == 0) { S = gen_0; P = gen_1; }
    1031             :   else
    1032             :   {
    1033             :     long j;
    1034         854 :     q = ginv(s); S = q; P = s;
    1035       16926 :     for (j = 1; j < k; j++)
    1036             :     {
    1037       16072 :       GEN sj = gaddgs(s, j);
    1038       16072 :       q = gmul(q, gdiv(x, sj));
    1039       16072 :       S = gadd(S, q);
    1040       16072 :       P = gmul(P, sj);
    1041             :     }
    1042         854 :     cox = expmx_xs(s, x, logx, prec); /* x^s exp(-x) */
    1043         854 :     S = gmul(S, gneg(cox));
    1044             :   }
    1045        2401 :   if (k && gequal0(sk))
    1046         175 :     return gadd(S, gdiv(eint1(x, prec), P));
    1047        2226 :   esk = gexpo(sk);
    1048        2226 :   if (esk > -7)
    1049             :   {
    1050        1015 :     GEN a, b, PG = gmul(sk, P);
    1051        1015 :     if (g) g = gmul(g, PG);
    1052        1015 :     a = incgam0(gaddgs(sk,1), x, g, prec);
    1053        1015 :     if (k == 0) cox = expmx_xs(s, x, logx, prec);
    1054        1015 :     b = gmul(gpowgs(x, k), cox);
    1055        1015 :     return gadd(S, gdiv(gsub(a, b), PG));
    1056             :   }
    1057        1211 :   E = prec2nbits(prec) + 1;
    1058        1211 :   if (gexpo(x) > 0)
    1059             :   {
    1060         420 :     long X = (long)(dblmodulus(x)/M_LN2);
    1061         420 :     prec += 2*nbits2extraprec(X);
    1062         420 :     x = gtofp(x, prec); mx = gneg(x);
    1063         420 :     logx = glog(x, prec); sk = gtofp(sk, prec);
    1064         420 :     E += X;
    1065             :   }
    1066        1211 :   if (isinexactreal(sk)) sk = gtofp(sk, prec+EXTRAPREC64);
    1067             :   /* |sk| < 2^-7 is small, guard against cancellation */
    1068        1211 :   F3 = gexpm1(gmul(sk, logx), prec);
    1069             :   /* ( gamma(1+sk) - exp(sk log(x))) ) / sk */
    1070        1211 :   S1 = gdiv(gsub(ggamma1m1(sk, prec+EXTRAPREC64), F3), sk);
    1071        1211 :   q = x; S3 = gdiv(x, gaddsg(1,sk));
    1072      255523 :   for (n = 2; gexpo(q) - gexpo(S3) > -E; n++)
    1073             :   {
    1074      254312 :     q = gmul(q, gdivgu(mx, n));
    1075      254312 :     S3 = gadd(S3, gdiv(q, gaddsg(n, sk)));
    1076             :   }
    1077        1211 :   S2 = gadd(gadd(S1, S3), gmul(F3, S3));
    1078        1211 :   return gadd(S, gdiv(S2, P));
    1079             : }
    1080             : 
    1081             : /* return |x| */
    1082             : double
    1083     9419000 : dblmodulus(GEN x)
    1084             : {
    1085     9419000 :   if (typ(x) == t_COMPLEX)
    1086             :   {
    1087     1505049 :     double a = gtodouble(gel(x,1));
    1088     1505049 :     double b = gtodouble(gel(x,2));
    1089     1505049 :     return sqrt(a*a + b*b);
    1090             :   }
    1091             :   else
    1092     7913951 :     return fabs(gtodouble(x));
    1093             : }
    1094             : 
    1095             : /* Driver routine. If g != NULL, assume that g=gamma(s,prec). */
    1096             : GEN
    1097       11543 : incgam0(GEN s, GEN x, GEN g, long prec)
    1098             : {
    1099             :   pari_sp av;
    1100             :   long E, l;
    1101             :   GEN z, rs, is;
    1102             : 
    1103       11543 :   if (gequal0(x)) return g? gcopy(g): ggamma(s,prec);
    1104       11543 :   if (gequal0(s)) return eint1(x, prec);
    1105        9737 :   l = precision(s); if (!l) l = prec;
    1106        9737 :   E = prec2nbits(l);
    1107        9737 :   if (gamma_use_asymp(x, E) ||
    1108        8316 :       (typ(s) == t_INT && signe(s) > 0 && gexpo(x) >= expi(s)))
    1109        2716 :     if ((z = incgam_asymp(s, x, l))) return z;
    1110        7140 :   av = avma; E++;
    1111        7140 :   rs = real_i(s);
    1112        7140 :   is = imag_i(s);
    1113             : #ifdef INCGAM_CF
    1114             :   /* Can one use continued fraction ? */
    1115             :   if (gequal0(is) && gequal0(imag_i(x)) && gsigne(x) > 0)
    1116             :   {
    1117             :     double sd = gtodouble(rs), LB, UB;
    1118             :     double xd = gtodouble(real_i(x));
    1119             :     if (sd > 0) {
    1120             :       LB = 15 + 0.1205*E;
    1121             :       UB = 5 + 0.752*E;
    1122             :     } else {
    1123             :       LB = -6 + 0.1205*E;
    1124             :       UB = 5 + 0.752*E + fabs(sd)/54.;
    1125             :     }
    1126             :     if (xd >= LB && xd <= UB)
    1127             :     {
    1128             :       if (DEBUGLEVEL > 2) err_printf("incgam: using continued fraction\n");
    1129             :       return gerepileupto(av, incgam_cf(s, x, xd, prec));
    1130             :     }
    1131             :   }
    1132             : #endif
    1133        7140 :   if (gsigne(rs) > 0 && gexpo(rs) >= -1)
    1134             :   { /* use complementary incomplete gamma */
    1135        4739 :     long n, egs, exd, precg, es = gexpo(s);
    1136        4739 :     if (es < 0) {
    1137         595 :       l += nbits2extraprec(-es) + 1;
    1138         595 :       x = gtofp(x, l);
    1139         595 :       if (isinexactreal(s)) s = gtofp(s, l);
    1140             :     }
    1141        4739 :     n = itos(gceil(rs));
    1142        4739 :     if (n > 100)
    1143             :     {
    1144             :       GEN gasx;
    1145         546 :       n -= 100;
    1146         546 :       if (es > 0)
    1147             :       {
    1148         546 :         es = mygamma(gtodouble(rs) - n, gtodouble(is)) / M_LN2;
    1149         546 :         if (es > 0)
    1150             :         {
    1151         546 :           l += nbits2extraprec(es);
    1152         546 :           x = gtofp(x, l);
    1153         546 :           if (isinexactreal(s)) s = gtofp(s, l);
    1154             :         }
    1155             :       }
    1156         546 :       gasx = incgam0(gsubgs(s, n), x, NULL, prec);
    1157         546 :       return gerepileupto(av, incgam_asymp_partial(s, x, gasx, n, prec));
    1158             :     }
    1159        4193 :     if (DEBUGLEVEL > 2) err_printf("incgam: using power series 1\n");
    1160             :     /* egs ~ expo(gamma(s)) */
    1161        4193 :     precg = g? precision(g): 0;
    1162        4193 :     egs = g? gexpo(g): (long)(mygamma(gtodouble(rs), gtodouble(is)) / M_LN2);
    1163        4193 :     if (egs > 0) {
    1164        1946 :       l += nbits2extraprec(egs) + 1;
    1165        1946 :       x = gtofp(x, l);
    1166        1946 :       if (isinexactreal(s)) s = gtofp(s, l);
    1167        1946 :       if (precg < l) g = NULL;
    1168             :     }
    1169        4193 :     z = incgamc_i(s, x, &exd, l);
    1170        4193 :     if (exd > 0)
    1171             :     {
    1172         896 :       l += nbits2extraprec(exd);
    1173         896 :       if (isinexactreal(s)) s = gtofp(s, l);
    1174         896 :       if (precg < l) g = NULL;
    1175             :     }
    1176             :     else
    1177             :     { /* gamma(s) negligible ? Compute to lower accuracy */
    1178        3297 :       long e = gexpo(z) - egs;
    1179        3297 :       if (e > 3)
    1180             :       {
    1181         420 :         E -= e;
    1182         420 :         if (E <= 0) g = gen_0; else if (!g) g = ggamma(s, nbits2prec(E));
    1183             :       }
    1184             :     }
    1185             :     /* worry about possible cancellation */
    1186        4193 :     if (!g) g = ggamma(s, maxss(l,precision(z)));
    1187        4193 :     return gerepileupto(av, gsub(g,z));
    1188             :   }
    1189        2401 :   if (DEBUGLEVEL > 2) err_printf("incgam: using power series 2\n");
    1190        2401 :   return gerepileupto(av, incgamspec(s, x, g, l));
    1191             : }
    1192             : 
    1193             : GEN
    1194        1106 : incgam(GEN s, GEN x, long prec) { return incgam0(s, x, NULL, prec); }
    1195             : 
    1196             : /* x a t_REAL */
    1197             : GEN
    1198        2156 : mpeint1(GEN x, GEN expx)
    1199             : {
    1200        2156 :   long s = signe(x);
    1201             :   pari_sp av;
    1202             :   GEN z;
    1203        2156 :   if (!s) pari_err_DOMAIN("eint1", "x","=",gen_0, x);
    1204        2149 :   if (s < 0) return eint1m(x, expx);
    1205        2009 :   z = cgetr(realprec(x));
    1206        2009 :   av = avma; affrr(eint1p(x, expx), z);
    1207        2009 :   set_avma(av); return z;
    1208             : }
    1209             : 
    1210             : static GEN
    1211         357 : cxeint1(GEN x, long prec)
    1212             : {
    1213         357 :   pari_sp av = avma, av2;
    1214             :   GEN q, S, run, z, H;
    1215         357 :   long n, E = prec2nbits(prec);
    1216             : 
    1217         357 :   if (gamma_use_asymp(x, E) && (z = eint1_asymp(x, NULL, prec))) return z;
    1218         252 :   E++;
    1219         252 :   if (gexpo(x) > 0)
    1220             :   { /* take cancellation into account, log2(\sum |x|^n / n!) = |x| / log(2) */
    1221          42 :     double dbx = dblmodulus(x);
    1222          42 :     long X = (long)((dbx + log(dbx))/M_LN2 + 10);
    1223          42 :     prec += nbits2extraprec(X);
    1224          42 :     x = gtofp(x, prec); E += X;
    1225             :   }
    1226         252 :   if (DEBUGLEVEL > 2) err_printf("eint1: using power series\n");
    1227         252 :   run = real_1(prec);
    1228         252 :   av2 = avma;
    1229         252 :   S = z = q = H = run;
    1230       48384 :   for (n = 2; gexpo(q) - gexpo(S) >= -E; n++)
    1231             :   {
    1232       48132 :     H = addrr(H, divru(run, n)); /* H = sum_{k<=n} 1/k */
    1233       48132 :     z = gdivgu(gmul(x,z), n);   /* z = x^(n-1)/n! */
    1234       48132 :     q = gmul(z, H); S = gadd(S, q);
    1235       48132 :     if ((n & 0x1ff) == 0) gerepileall(av2, 4, &z, &q, &S, &H);
    1236             :   }
    1237         252 :   S = gmul(gmul(x, S), gexp(gneg_i(x), prec));
    1238         252 :   return gerepileupto(av, gsub(S, gadd(glog(x, prec), mpeuler(prec))));
    1239             : }
    1240             : 
    1241             : GEN
    1242        2513 : eint1(GEN x, long prec)
    1243             : {
    1244        2513 :   switch(typ(x))
    1245             :   {
    1246         357 :     case t_COMPLEX: return cxeint1(x, prec);
    1247        1771 :     case t_REAL: break;
    1248         385 :     default: x = gtofp(x, prec);
    1249             :   }
    1250        2156 :   return mpeint1(x,NULL);
    1251             : }
    1252             : 
    1253             : GEN
    1254          49 : veceint1(GEN C, GEN nmax, long prec)
    1255             : {
    1256          49 :   if (!nmax) return eint1(C,prec);
    1257           7 :   if (typ(nmax) != t_INT) pari_err_TYPE("veceint1",nmax);
    1258           7 :   if (typ(C) != t_REAL) {
    1259           7 :     C = gtofp(C, prec);
    1260           7 :     if (typ(C) != t_REAL) pari_err_TYPE("veceint1",C);
    1261             :   }
    1262           7 :   if (signe(C) <= 0) pari_err_DOMAIN("veceint1", "argument", "<=", gen_0,C);
    1263           7 :   return mpveceint1(C, NULL, itos(nmax));
    1264             : }
    1265             : 
    1266             : /* j > 0, a t_REAL. Return sum_{m >= 0} a^m / j(j+1)...(j+m)).
    1267             :  * Stop when expo(summand) < E; note that s(j-1) = (a s(j) + 1) / (j-1). */
    1268             : static GEN
    1269         259 : mp_sum_j(GEN a, long j, long E, long prec)
    1270             : {
    1271         259 :   pari_sp av = avma;
    1272         259 :   GEN q = divru(real_1(prec), j), s = q;
    1273             :   long m;
    1274         259 :   for (m = 0;; m++)
    1275             :   {
    1276        4871 :     if (expo(q) < E) break;
    1277        4612 :     q = mulrr(q, divru(a, m+j));
    1278        4612 :     s = addrr(s, q);
    1279             :   }
    1280         259 :   return gerepileuptoleaf(av, s);
    1281             : }
    1282             : /* Return the s_a(j), j <= J */
    1283             : static GEN
    1284         259 : sum_jall(GEN a, long J, long prec)
    1285             : {
    1286         259 :   GEN s = cgetg(J+1, t_VEC);
    1287         259 :   long j, E = -prec2nbits(prec) - 5;
    1288         259 :   gel(s, J) = mp_sum_j(a, J, E, prec);
    1289       11103 :   for (j = J-1; j; j--)
    1290       10844 :     gel(s,j) = divru(addrs(mulrr(a, gel(s,j+1)), 1), j);
    1291         259 :   return s;
    1292             : }
    1293             : 
    1294             : /* T a dense t_POL with t_REAL coeffs. Return T(n) [faster than poleval] */
    1295             : static GEN
    1296      481831 : rX_s_eval(GEN T, long n)
    1297             : {
    1298      481831 :   long i = lg(T)-1;
    1299      481831 :   GEN c = gel(T,i);
    1300    10743228 :   for (i--; i>=2; i--) c = gadd(mulrs(c,n),gel(T,i));
    1301      481831 :   return c;
    1302             : }
    1303             : 
    1304             : /* C>0 t_REAL, eC = exp(C). Return eint1(n*C) for 1<=n<=N. Absolute accuracy */
    1305             : GEN
    1306         266 : mpveceint1(GEN C, GEN eC, long N)
    1307             : {
    1308         266 :   const long prec = realprec(C);
    1309         266 :   long Nmin = 15; /* >= 1. E.g. between 10 and 30, but little effect */
    1310         266 :   GEN en, v, w = cgetg(N+1, t_VEC);
    1311             :   pari_sp av0;
    1312             :   double DL;
    1313             :   long n, j, jmax, jmin;
    1314         266 :   if (!N) return w;
    1315      486017 :   for (n = 1; n <= N; n++) gel(w,n) = cgetr(prec);
    1316         266 :   av0 = avma;
    1317         266 :   if (N < Nmin) Nmin = N;
    1318         266 :   if (!eC) eC = mpexp(C);
    1319         266 :   en = eC; affrr(eint1p(C, en), gel(w,1));
    1320        3920 :   for (n = 2; n <= Nmin; n++)
    1321             :   {
    1322             :     pari_sp av2;
    1323        3654 :     en = mulrr(en,eC); /* exp(n C) */
    1324        3654 :     av2 = avma;
    1325        3654 :     affrr(eint1p(mulru(C,n), en), gel(w,n));
    1326        3654 :     set_avma(av2);
    1327             :   }
    1328         266 :   if (Nmin == N) { set_avma(av0); return w; }
    1329             : 
    1330         259 :   DL = prec2nbits_mul(prec, M_LN2) + 5;
    1331         259 :   jmin = ceil(DL/log((double)N)) + 1;
    1332         259 :   jmax = ceil(DL/log((double)Nmin)) + 1;
    1333         259 :   v = sum_jall(C, jmax, prec);
    1334         259 :   en = powrs(eC, -N); /* exp(-N C) */
    1335         259 :   affrr(eint1p(mulru(C,N), invr(en)), gel(w,N));
    1336        7038 :   for (j = jmin, n = N-1; j <= jmax; j++)
    1337             :   {
    1338        6779 :     long limN = maxss((long)ceil(exp(DL/j)), Nmin);
    1339             :     GEN polsh;
    1340        6779 :     setlg(v, j+1);
    1341        6779 :     polsh = RgV_to_RgX_reverse(v, 0);
    1342      488610 :     for (; n >= limN; n--)
    1343             :     {
    1344      481831 :       pari_sp av2 = avma;
    1345      481831 :       GEN S = divri(mulrr(en, rX_s_eval(polsh, -n)), powuu(n,j));
    1346             :       /* w[n+1] - exp(-n C) * polsh(-n) / (-n)^j */
    1347      481831 :       GEN c = odd(j)? addrr(gel(w,n+1), S) : subrr(gel(w,n+1), S);
    1348      481831 :       affrr(c, gel(w,n)); set_avma(av2);
    1349      481831 :       en = mulrr(en,eC); /* exp(-n C) */
    1350             :     }
    1351             :   }
    1352         259 :   set_avma(av0); return w;
    1353             : }
    1354             : 
    1355             : /* erfc via numerical integration : assume real(x)>=1 */
    1356             : static GEN
    1357          14 : cxerfc_r1(GEN x, long prec)
    1358             : {
    1359             :   GEN h, h2, eh2, denom, res, lambda;
    1360             :   long u, v;
    1361          14 :   const double D = prec2nbits_mul(prec, M_LN2);
    1362          14 :   const long npoints = (long)ceil(D/M_PI)+1;
    1363          14 :   pari_sp av = avma;
    1364             :   {
    1365          14 :     double t = exp(-2*M_PI*M_PI/D); /* ~exp(-2*h^2) */
    1366          14 :     v = 30; /* bits that fit in both long and double mantissa */
    1367          14 :     u = (long)floor(t*(1L<<v));
    1368             :     /* define exp(-2*h^2) to be u*2^(-v) */
    1369             :   }
    1370          14 :   incrprec(prec);
    1371          14 :   x = gtofp(x,prec);
    1372          14 :   eh2 = sqrtr_abs(rtor(shiftr(dbltor(u),-v),prec));
    1373          14 :   h2 = negr(logr_abs(eh2));
    1374          14 :   h = sqrtr_abs(h2);
    1375          14 :   lambda = gdiv(x,h);
    1376          14 :   denom = gsqr(lambda);
    1377             :   { /* res = h/x + 2*x*h*sum(k=1,npoints,exp(-(k*h)^2)/(lambda^2+k^2)); */
    1378             :     GEN Uk; /* = exp(-(kh)^2) */
    1379          14 :     GEN Vk = eh2;/* = exp(-(2k+1)h^2) */
    1380          14 :     pari_sp av2 = avma;
    1381             :     long k;
    1382             :     /* k = 0 moved out for efficiency */
    1383          14 :     denom = gaddsg(1,denom);
    1384          14 :     Uk = Vk;
    1385          14 :     Vk = mulur(u,Vk); shiftr_inplace(Vk, -v);
    1386          14 :     res = gdiv(Uk, denom);
    1387         420 :     for (k = 1; k < npoints; k++)
    1388             :     {
    1389         406 :       if ((k & 255) == 0) gerepileall(av2,4,&denom,&Uk,&Vk,&res);
    1390         406 :       denom = gaddsg(2*k+1,denom);
    1391         406 :       Uk = mpmul(Uk,Vk);
    1392         406 :       Vk = mulur(u,Vk); shiftr_inplace(Vk, -v);
    1393         406 :       res = gadd(res, gdiv(Uk, denom));
    1394             :     }
    1395             :   }
    1396          14 :   res = gmul(res, gshift(lambda,1));
    1397             :   /* 0 term : */
    1398          14 :   res = gadd(res, ginv(lambda));
    1399          14 :   res = gmul(res, gdiv(gexp(gneg(gsqr(x)), prec), mppi(prec)));
    1400          14 :   if (rtodbl(real_i(x)) < sqrt(D))
    1401             :   {
    1402          14 :     GEN t = gmul(divrr(Pi2n(1,prec),h), x);
    1403          14 :     res = gsub(res, gdivsg(2, cxexpm1(t, prec)));
    1404             :   }
    1405          14 :   return gerepileupto(av,res);
    1406             : }
    1407             : 
    1408             : GEN
    1409          56 : gerfc(GEN x, long prec)
    1410             : {
    1411             :   GEN z, xr, xi, res;
    1412             :   long s;
    1413             :   pari_sp av;
    1414             : 
    1415          56 :   x = trans_fix_arg(&prec,&x,&xr,&xi, &av,&res);
    1416          56 :   s = signe(xr);
    1417          56 :   if (s > 0 || (s == 0 && signe(xi) >= 0)) {
    1418          84 :     if (cmprs(xr, 1) > 0) /* use numerical integration */
    1419          14 :       z = cxerfc_r1(x, prec);
    1420             :     else
    1421             :     { /* erfc(x) = incgam(1/2,x^2)/sqrt(Pi) */
    1422          28 :       GEN sqrtpi = sqrtr(mppi(prec));
    1423          28 :       z = incgam0(ghalf, gsqr(x), sqrtpi, prec);
    1424          28 :       z = gdiv(z, sqrtpi);
    1425             :     }
    1426             :   }
    1427             :   else
    1428             :   { /* erfc(-x)=2-erfc(x) */
    1429             :     /* FIXME could decrease prec
    1430             :     long size = nbits2extraprec((imag(x)^2-real(x)^2)/log(2));
    1431             :     prec = size > 0 ? prec : prec + size;
    1432             :     */
    1433             :     /* NOT gsubsg(2, ...) : would create a result of
    1434             :      * huge accuracy if re(x)>>1, rounded to 2 by subsequent affc_fixlg... */
    1435          14 :     z = gsub(real2n(1,prec+EXTRAPREC64), gerfc(gneg(x), prec));
    1436             :   }
    1437          56 :   set_avma(av); return affc_fixlg(z, res);
    1438             : }
    1439             : 
    1440             : /***********************************************************************/
    1441             : /**                                                                   **/
    1442             : /**                      RIEMANN ZETA FUNCTION                        **/
    1443             : /**                                                                   **/
    1444             : /***********************************************************************/
    1445             : static const double log2PI = 1.83787706641;
    1446             : 
    1447             : static double
    1448        4578 : get_xinf(double beta)
    1449             : {
    1450        4578 :   const double maxbeta = 0.06415003; /* 3^(-2.5) */
    1451             :   double x0, y0, x1;
    1452             : 
    1453        4578 :   if (beta < maxbeta) return beta + pow(3*beta, 1.0/3.0);
    1454        4578 :   x0 = beta + M_PI/2.0;
    1455             :   for(;;)
    1456             :   {
    1457       10486 :     y0 = x0*x0;
    1458        7532 :     x1 = (beta+atan(x0)) * (1+y0) / y0 - 1/x0;
    1459        7532 :     if (0.99*x0 < x1) return x1;
    1460        2954 :     x0 = x1;
    1461             :   }
    1462             : }
    1463             : /* optimize for zeta( s + it, prec ), assume |s-1| > 0.1
    1464             :  * (if gexpo(u = s-1) < -5, we use the functional equation s->1-s) */
    1465             : static void
    1466        8575 : optim_zeta(GEN S, long prec, long *pp, long *pn)
    1467             : {
    1468             :   double s, t, alpha, beta, n, B;
    1469             :   long p;
    1470        8575 :   if (typ(S) == t_REAL) {
    1471         266 :     s = rtodbl(S);
    1472         266 :     t = 0.;
    1473             :   } else {
    1474        8309 :     s = rtodbl(gel(S,1));
    1475        8309 :     t = fabs( rtodbl(gel(S,2)) );
    1476             :   }
    1477             : 
    1478        8575 :   B = prec2nbits_mul(prec, M_LN2);
    1479        8575 :   if (s > 0 && !t) /* positive real input */
    1480             :   {
    1481         238 :     beta = B + 0.61 + s*(log2PI - log(s));
    1482         476 :     if (beta > 0)
    1483             :     {
    1484         231 :       p = (long)ceil(beta / 2.0);
    1485         231 :       n = fabs(s + 2*p-1)/(2*M_PI);
    1486             :     }
    1487             :     else
    1488             :     {
    1489           7 :       p = 0;
    1490           7 :       n = exp((B - M_LN2) / s);
    1491             :     }
    1492             :   }
    1493        8337 :   else if (s <= 0 || t < 0.01) /* s < 0 may occur if s ~ 0 */
    1494        3752 :   { /* TODO: the crude bounds below are generally valid. Optimize ? */
    1495        3752 :     double l,l2, la = 1.; /* heuristic */
    1496        3752 :     double rlog, ilog; dcxlog(s-1,t, &rlog,&ilog);
    1497        3752 :     l2 = (s - 0.5)*rlog - t*ilog; /* = Re( (S - 1/2) log (S-1) ) */
    1498        3752 :     l = (B - l2 + s*log2PI) / (2. * (1.+ log((double)la)));
    1499        3752 :     l2 = dabs(s, t)/2;
    1500        3752 :     if (l < l2) l = l2;
    1501        3752 :     p = (long) ceil(l); if (p < 2) p = 2;
    1502        3752 :     n = 1 + dabs(p+s/2.-.25, t/2) * la / M_PI;
    1503             :   }
    1504             :   else
    1505             :   {
    1506        4585 :     double sn = dabs(s, t), L = log(sn/s);
    1507        4585 :     alpha = B - 0.39 + L + s*(log2PI - log(sn));
    1508        4585 :     beta = (alpha+s)/t - atan(s/t);
    1509        4585 :     p = 0;
    1510        4585 :     if (beta > 0)
    1511             :     {
    1512        4578 :       beta = 1.0 - s + t * get_xinf(beta);
    1513        4578 :       if (beta > 0) p = (long)ceil(beta / 2.0);
    1514             :     }
    1515             :     else
    1516           7 :       if (s < 1.0) p = 1;
    1517        4585 :     n = p? dabs(s + 2*p-1, t) / (2*M_PI) : exp((B-M_LN2+L) / s);
    1518             :   }
    1519        8575 :   *pp = p;
    1520        8575 :   *pn = (long)ceil(n);
    1521        8575 :   if (*pp < 0 || *pn < 0) pari_err_OVERFLOW("zeta");
    1522        8575 : }
    1523             : 
    1524             : /* zeta(a*j+b), j=0..N-1, b>1, using sumalt. Johansonn's thesis, Algo 4.7.1 */
    1525             : static GEN
    1526         723 : veczetas(long a, long b, long N, long prec)
    1527             : {
    1528         723 :   const long n = ceil(2 + prec2nbits_mul(prec, M_LN2/1.7627));
    1529         723 :   pari_sp av = avma;
    1530         723 :   GEN c, d, z = zerovec(N);
    1531             :   long j, k;
    1532         723 :   c = d = int2n(2*n-1);
    1533       85914 :   for (k = n; k > 1; k--)
    1534             :   {
    1535       85191 :     GEN u, t = divii(d, powuu(k,b));
    1536       85193 :     if (!odd(k)) t = negi(t);
    1537       85347 :     gel(z,1) = addii(gel(z,1), t);
    1538       85105 :     u = powuu(k,a);
    1539     3917942 :     for (j = 1; j < N; j++)
    1540             :     {
    1541     3877840 :       t = divii(t,u); if (!signe(t)) break;
    1542     3835536 :       gel(z,j+1) = addii(gel(z,j+1), t);
    1543             :     }
    1544       85085 :     c = muluui(k,2*k-1,c);
    1545       85102 :     c = diviuuexact(c, 2*(n-k+1),n+k-1);
    1546       85146 :     d = addii(d,c);
    1547       85130 :     if (gc_needed(av,3))
    1548             :     {
    1549           9 :       if(DEBUGMEM>1) pari_warn(warnmem,"zetaBorwein, k = %ld", k);
    1550           9 :       gerepileall(av, 3, &c,&d,&z);
    1551             :     }
    1552             :   }
    1553             :   /* k = 1 */
    1554       52436 :   for (j = 1; j <= N; j++) gel(z,j) = addii(gel(z,j), d);
    1555         692 :   d = addiu(d, 1);
    1556       52687 :   for (j = 0, k = b - 1; j < N; j++, k += a)
    1557       51964 :     gel(z,j+1) = rdivii(shifti(gel(z,j+1), k), subii(shifti(d,k), d), prec);
    1558         723 :   return z;
    1559             : }
    1560             : /* zeta(a*j+b), j=0..N-1, b>1, using sumalt */
    1561             : GEN
    1562         217 : veczeta(GEN a, GEN b, long N, long prec)
    1563             : {
    1564         217 :   pari_sp av = avma;
    1565             :   long n, j, k;
    1566         217 :   GEN L, c, d, z = zerovec(N);
    1567         217 :   if (typ(a) == t_INT && typ(b) == t_INT)
    1568         126 :     return gerepilecopy(av, veczetas(itos(a),  itos(b), N, prec));
    1569          91 :   n = ceil(2 + prec2nbits_mul(prec, M_LN2/1.7627));
    1570          91 :   c = d = int2n(2*n-1);
    1571       13819 :   for (k = n; k; k--)
    1572             :   {
    1573             :     GEN u, t;
    1574       13728 :     L = logr_abs(utor(k, prec)); /* log(k) */
    1575       13728 :     t = gdiv(d, gexp(gmul(b, L), prec)); /* d / k^b */
    1576       13728 :     if (!odd(k)) t = gneg(t);
    1577       13728 :     gel(z,1) = gadd(gel(z,1), t);
    1578       13728 :     u = gexp(gmul(a, L), prec);
    1579      896888 :     for (j = 1; j < N; j++)
    1580             :     {
    1581      889583 :       t = gdiv(t,u); if (gexpo(t) < 0) break;
    1582      883160 :       gel(z,j+1) = gadd(gel(z,j+1), t);
    1583             :     }
    1584       13728 :     c = muluui(k,2*k-1,c);
    1585       13728 :     c = diviuuexact(c, 2*(n-k+1),n+k-1);
    1586       13728 :     d = addii(d,c);
    1587       13728 :     if (gc_needed(av,3))
    1588             :     {
    1589           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"veczeta, k = %ld", k);
    1590           0 :       gerepileall(av, 3, &c,&d,&z);
    1591             :     }
    1592             :   }
    1593          91 :   L = mplog2(prec);
    1594        5796 :   for (j = 0; j < N; j++)
    1595             :   {
    1596        5705 :     GEN u = gsubgs(gadd(b, gmulgu(a,j)), 1);
    1597        5705 :     GEN w = gexp(gmul(u, L), prec); /* 2^u */
    1598        5705 :     gel(z,j+1) = gdiv(gmul(gel(z,j+1), w), gmul(d,gsubgs(w,1)));
    1599             :   }
    1600          91 :   return gerepilecopy(av, z);
    1601             : }
    1602             : 
    1603             : GEN
    1604       26310 : constzeta(long n, long prec)
    1605             : {
    1606       26310 :   GEN o = zetazone, z;
    1607       26310 :   long l = o? lg(o): 0;
    1608             :   pari_sp av;
    1609       26310 :   if (l > n)
    1610             :   {
    1611       25864 :     long p = realprec(gel(o,1));
    1612       25864 :     if (p >= prec) return o;
    1613             :   }
    1614         597 :   n = maxss(n, l + 15);
    1615         597 :   av = avma; z = veczetas(1, 2, n-1, prec);
    1616         597 :   zetazone = gclone(vec_prepend(z, mpeuler(prec)));
    1617         597 :   set_avma(av); guncloneNULL(o); return zetazone;
    1618             : }
    1619             : 
    1620             : /* zeta(s) using sumalt, case h=0,N=1. Assume s > 1 */
    1621             : static GEN
    1622         598 : zetaBorwein(long s, long prec)
    1623             : {
    1624         598 :   pari_sp av = avma;
    1625         598 :   const long n = ceil(2 + prec2nbits_mul(prec, M_LN2/1.7627));
    1626             :   long k;
    1627         598 :   GEN c, d, z = gen_0;
    1628         598 :   c = d = int2n(2*n-1);
    1629      135672 :   for (k = n; k; k--)
    1630             :   {
    1631      135074 :     GEN t = divii(d, powuu(k,s));
    1632      135074 :     z = odd(k)? addii(z,t): subii(z,t);
    1633      135074 :     c = muluui(k,2*k-1,c);
    1634      135074 :     c = diviuuexact(c, 2*(n-k+1),n+k-1);
    1635      135074 :     d = addii(d,c);
    1636      135074 :     if (gc_needed(av,3))
    1637             :     {
    1638           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"zetaBorwein, k = %ld", k);
    1639           0 :       gerepileall(av, 3, &c,&d,&z);
    1640             :     }
    1641             :   }
    1642         598 :   return rdivii(shifti(z, s-1), subii(shifti(d,s-1), d), prec);
    1643             : }
    1644             : 
    1645             : /* assume k != 1 */
    1646             : GEN
    1647        4739 : szeta(long k, long prec)
    1648             : {
    1649        4739 :   pari_sp av = avma;
    1650             :   GEN z;
    1651             : 
    1652        4739 :   if (!k) { z = real2n(-1, prec); setsigne(z,-1); return z; }
    1653        4732 :   if (k < 0)
    1654             :   {
    1655         252 :     if (!odd(k)) return gen_0;
    1656             :     /* the one value such that k < 0 and 1 - k < 0, due to overflow */
    1657         252 :     if ((ulong)k == (HIGHBIT | 1))
    1658           0 :       pari_err_OVERFLOW("zeta [large negative argument]");
    1659         252 :     k = 1-k;
    1660         252 :     z = bernreal(k, prec); togglesign(z);
    1661         252 :     return gerepileuptoleaf(av, divru(z, k));
    1662             :   }
    1663             :   /* k > 1 */
    1664        4480 :   if (k > prec2nbits(prec)+1) return real_1(prec);
    1665        4473 :   if (zetazone && realprec(gel(zetazone,1)) >= prec && lg(zetazone) > k)
    1666         119 :     return rtor(gel(zetazone, k), prec);
    1667        4354 :   if (!odd(k))
    1668             :   {
    1669             :     GEN B;
    1670        2583 :     if (!bernzone) constbern(0);
    1671        2583 :     if (k < lg(bernzone))
    1672        2128 :       B = gel(bernzone, k>>1);
    1673             :     else
    1674             :     {
    1675         455 :       if (bernbitprec(k) > prec2nbits(prec))
    1676           0 :         return gerepileupto(av, invr(inv_szeta_euler(k, prec)));
    1677         455 :       B = bernfrac(k);
    1678             :     }
    1679             :     /* B = B_k */
    1680        2583 :     z = gmul(powru(Pi2n(1, prec + EXTRAPRECWORD), k), B);
    1681        2583 :     z = divrr(z, mpfactr(k,prec));
    1682        2583 :     setsigne(z, 1); shiftr_inplace(z, -1);
    1683             :   }
    1684             :   else
    1685             :   {
    1686        1771 :     double p = prec2nbits_mul(prec,0.393); /* bit / log_2(3+sqrt(8)) */
    1687        1771 :     p = log2(p * log(p));
    1688        3542 :     z = (p * k > prec2nbits(prec))? invr(inv_szeta_euler(k, prec))
    1689        1771 :                                   : zetaBorwein(k, prec);
    1690             :   }
    1691        4354 :   return gerepileuptoleaf(av, z);
    1692             : }
    1693             : 
    1694             : /* Ensure |1-s| >= 1/32 and (|s| <= 1/32 or real(s) >= 1/2) */
    1695             : static int
    1696        8589 : zeta_funeq(GEN *ps)
    1697             : {
    1698        8589 :   GEN s = *ps, u;
    1699        8589 :   if (typ(s) == t_REAL)
    1700             :   {
    1701         273 :     u = subsr(1, s);
    1702         273 :     if (expo(u) >= -5
    1703         266 :         && ((signe(s) > 0 && expo(s) >= -1) || expo(s) <= -5)) return 0;
    1704             :   }
    1705             :   else
    1706             :   {
    1707        8316 :     GEN sig = gel(s,1);
    1708        8316 :     u = gsubsg(1, s);
    1709        8316 :     if (gexpo(u) >= -5
    1710        8309 :         && ((signe(sig) > 0 && expo(sig) >= -1) || gexpo(s) <= -5)) return 0;
    1711             :   }
    1712        1505 :   *ps = u; return 1;
    1713             : }
    1714             : /* s0 a t_INT, t_REAL or t_COMPLEX.
    1715             :  * If a t_INT, assume it's not a trivial case (i.e we have s0 > 1, odd) */
    1716             : static GEN
    1717        8596 : czeta(GEN s0, long prec)
    1718             : {
    1719        8596 :   GEN ms, s, u, y, res, tes, sig, tau, invn2, ns, Ns, funeq_factor = NULL;
    1720             :   long i, nn, lim, lim2;
    1721        8596 :   pari_sp av0 = avma, av, av2;
    1722             :   pari_timer T;
    1723             : 
    1724        8596 :   if (DEBUGLEVEL>2) timer_start(&T);
    1725        8596 :   s = trans_fix_arg(&prec,&s0,&sig,&tau,&av,&res);
    1726        8596 :   if (typ(s0) == t_INT) return gerepileupto(av0, gzeta(s0, prec));
    1727        8589 :   if (zeta_funeq(&s)) /* s -> 1-s */
    1728             :   { /* Gamma(s) (2Pi)^-s 2 cos(Pi s/2) [new s] */
    1729        1505 :     GEN t = gmul(ggamma(s,prec), pow2Pis(gsubgs(s0,1), prec));
    1730        1505 :     sig = real_i(s);
    1731        1505 :     funeq_factor = gmul2n(gmul(t, gsin(gmul(Pi2n(-1,prec),s0), prec)), 1);
    1732             :   }
    1733        8589 :   if (gcmpgs(sig, prec2nbits(prec) + 1) > 0) { /* zeta(s) = 1 */
    1734          14 :     if (!funeq_factor) { set_avma(av0); return real_1(prec); }
    1735           7 :     return gerepileupto(av0, funeq_factor);
    1736             :   }
    1737        8575 :   optim_zeta(s, prec, &lim, &nn);
    1738        8575 :   if (DEBUGLEVEL>2) err_printf("lim, nn: [%ld, %ld]\n", lim, nn);
    1739             : 
    1740        8575 :   ms = gneg(s);
    1741        8575 :   if (umuluu_le(nn, prec, 10000000))
    1742             :   {
    1743        8575 :     incrprec(prec); /* one extra word of precision */
    1744        8575 :     Ns = vecpowug(nn, ms, prec);
    1745        8575 :     ns = gel(Ns,nn); setlg(Ns, nn);
    1746        8575 :     y = gadd(gmul2n(ns, -1), RgV_sum(Ns));
    1747             :   }
    1748             :   else
    1749             :   {
    1750           0 :     Ns = dirpowerssum(nn, ms, prec);
    1751           0 :     incrprec(prec); /* one extra word of precision */
    1752           0 :     ns = gpow(utor(nn, prec), ms, prec);
    1753           0 :     y = gsub(Ns, gmul2n(ns, -1));
    1754             :   }
    1755        8575 :   if (DEBUGLEVEL>2) timer_printf(&T,"sum from 1 to N");
    1756        8575 :   constbern(lim);
    1757        8575 :   if (DEBUGLEVEL>2) timer_start(&T);
    1758        8575 :   invn2 = divri(real_1(prec), sqru(nn)); lim2 = lim<<1;
    1759        8575 :   tes = bernfrac(lim2);
    1760             :   {
    1761             :     GEN s1, s2, s3, s4, s5;
    1762        8575 :     s2 = gmul(s, gsubgs(s,1));
    1763        8575 :     s3 = gmul2n(invn2,3);
    1764        8575 :     av2 = avma;
    1765        8575 :     s1 = gsubgs(gmul2n(s,1), 1);
    1766        8575 :     s4 = gmul(invn2, gmul2n(gaddsg(4*lim-2,s1),1));
    1767        8575 :     s5 = gmul(invn2, gadd(s2, gmulsg(lim2, gaddgs(s1, lim2))));
    1768      805263 :     for (i = lim2-2; i>=2; i -= 2)
    1769             :     {
    1770      796688 :       s5 = gsub(s5, s4);
    1771      796688 :       s4 = gsub(s4, s3);
    1772      796688 :       tes = gadd(bernfrac(i), gdivgunextu(gmul(s5,tes), i+1));
    1773      796688 :       if (gc_needed(av2,3))
    1774             :       {
    1775           0 :         if(DEBUGMEM>1) pari_warn(warnmem,"czeta i = %ld", i);
    1776           0 :         gerepileall(av2,3, &tes,&s5,&s4);
    1777             :       }
    1778             :     }
    1779        8575 :     u = gmul(gmul(tes,invn2), gmul2n(s2, -1));
    1780        8575 :     tes = gmulsg(nn, gaddsg(1, u));
    1781             :   }
    1782        8575 :   if (DEBUGLEVEL>2) timer_printf(&T,"Bernoulli sum");
    1783             :   /* y += tes n^(-s) / (s-1) */
    1784        8575 :   y = gadd(y, gmul(tes, gdiv(ns, gsubgs(s,1))));
    1785        8575 :   if (funeq_factor) y = gmul(y, funeq_factor);
    1786        8575 :   set_avma(av); return affc_fixlg(y,res);
    1787             : }
    1788             : /* v a t_VEC/t_COL */
    1789             : int
    1790          14 : RgV_is_arithprog(GEN v, GEN *a, GEN *b)
    1791             : {
    1792          14 :   pari_sp av = avma, av2;
    1793          14 :   long i, n = lg(v)-1;
    1794          14 :   if (n == 0) { *a = *b = gen_0; return 1; }
    1795          14 :   *a = gel(v,1);
    1796          14 :   if (n == 1) { * b = gen_0; return 1; }
    1797          14 :   *b = gsub(gel(v,2), *a); av2 = avma;
    1798          28 :   for (i = 2; i < n; i++)
    1799          14 :     if (!gequal(*b, gsub(gel(v,i+1), gel(v,i)))) return gc_int(av,0);
    1800          14 :   return gc_int(av2,1);
    1801             : }
    1802             : 
    1803             : GEN
    1804       13349 : gzeta(GEN x, long prec)
    1805             : {
    1806       13349 :   pari_sp av = avma;
    1807             :   GEN y;
    1808       13349 :   if (gequal1(x)) pari_err_DOMAIN("zeta", "argument", "=", gen_1, x);
    1809       13328 :   switch(typ(x))
    1810             :   {
    1811        4501 :     case t_INT:
    1812        4501 :       if (is_bigint(x))
    1813             :       {
    1814          21 :         if (signe(x) > 0) return real_1(prec);
    1815          14 :         if (mod2(x) == 0) return real_0(prec);
    1816           7 :         pari_err_OVERFLOW("zeta [large negative argument]");
    1817             :       }
    1818        4480 :       return szeta(itos(x),prec);
    1819        8596 :     case t_REAL: case t_COMPLEX: return czeta(x,prec);
    1820          14 :     case t_PADIC: return Qp_zeta(x);
    1821          21 :     case t_VEC: case t_COL:
    1822             :     {
    1823             :       GEN a, b;
    1824          21 :       long n = lg(x) - 1;
    1825          21 :       if (n > 1 && RgV_is_arithprog(x, &a, &b))
    1826             :       {
    1827          14 :         if (!is_real_t(typ(a)) || !is_real_t(typ(b)) || gcmpgs(a, 1) <= 0)
    1828           0 :         { set_avma(av); break; }
    1829          14 :         a = veczeta(b, a, n, prec);
    1830          14 :         settyp(a, typ(x)); return a;
    1831             :       }
    1832             :     }
    1833             :     default:
    1834         203 :       if (!(y = toser_i(x))) break;
    1835          35 :       if (gequal1(y))
    1836           7 :         pari_err_DOMAIN("zeta", "argument", "=", gen_1, y);
    1837          28 :       return gerepileupto(av, lfun(gen_1,y,prec2nbits(prec)));
    1838             :   }
    1839         168 :   return trans_eval("zeta",gzeta,x,prec);
    1840             : }
    1841             : 
    1842             : /***********************************************************************/
    1843             : /**                                                                   **/
    1844             : /**                    FONCTIONS POLYLOGARITHME                       **/
    1845             : /**                                                                   **/
    1846             : /***********************************************************************/
    1847             : 
    1848             : /* smallish k such that bernbitprec(K) > bit + Kdz, K = 2k+4 */
    1849             : static long
    1850          21 : get_k(double dz, long bit)
    1851             : {
    1852             :   long a, b;
    1853          21 :   for (b = 128;; b <<= 1)
    1854          21 :     if (bernbitprec(b) > bit + b*dz) break;
    1855          21 :   if (b == 128) return 128;
    1856           0 :   a = b >> 1;
    1857           0 :   while (b - a > 64)
    1858             :   {
    1859           0 :     long c = (a+b) >> 1;
    1860           0 :     if (bernbitprec(c) > bit + c*dz) b = c; else a = c;
    1861             :   }
    1862           0 :   return b >> 1;
    1863             : }
    1864             : 
    1865             : /* m >= 2. Validity domain |log x| < 2*Pi, contains log |x| < 5.44,
    1866             :  * Li_m(x = e^z) = sum_{n >= 0} zeta(m-n) z^n / n!
    1867             :  *    with zeta(1) := H_{m-1} - log(-z) */
    1868             : static GEN
    1869          21 : cxpolylog(long m, GEN x, long prec)
    1870             : {
    1871             :   long li, n, k, ksmall, real;
    1872             :   GEN vz, z, Z, h, q, s, S;
    1873             :   pari_sp av;
    1874             :   double dz;
    1875             :   pari_timer T;
    1876             : 
    1877          21 :   if (gequal1(x)) return szeta(m,prec);
    1878             :   /* x real <= 1 ==> Li_m(x) real */
    1879          21 :   real = (typ(x) == t_REAL && (expo(x) < 0 || signe(x) <= 0));
    1880             : 
    1881          21 :   vz = constzeta(m, prec);
    1882          21 :   z = glog(x,prec);
    1883             :   /* n = 0 */
    1884          21 :   q = gen_1; s = gel(vz, m);
    1885          28 :   for (n=1; n < m-1; n++)
    1886             :   {
    1887           7 :     q = gdivgu(gmul(q,z),n);
    1888           7 :     s = gadd(s, gmul(gel(vz,m-n), real? real_i(q): q));
    1889             :   }
    1890             :   /* n = m-1 */
    1891          21 :     q = gdivgu(gmul(q,z),n); /* multiply by "zeta(1)" */
    1892          21 :     h = gmul(q, gsub(harmonic(m-1), glog(gneg_i(z),prec)));
    1893          21 :     s = gadd(s, real? real_i(h): h);
    1894             :   /* n = m */
    1895          21 :     q = gdivgu(gmul(q,z),m);
    1896          21 :     s = gadd(s, gdivgs(real? real_i(q): q, -2)); /* zeta(0) = -1/2 */
    1897             :   /* n = m+1 */
    1898          21 :     q = gdivgu(gmul(q,z),m+1); /* = z^(m+1) / (m+1)! */
    1899          21 :     s = gadd(s, gdivgs(real? real_i(q): q, -12)); /* zeta(-1) = -1/12 */
    1900             : 
    1901          21 :   li = -(prec2nbits(prec)+1);
    1902          21 :   if (DEBUGLEVEL) timer_start(&T);
    1903          21 :   dz = dbllog2(z) - log2PI; /*  ~ log2(|z|/2Pi) */
    1904             :   /* sum_{k >= 1} zeta(-1-2k) * z^(2k+m+1) / (2k+m+1)!
    1905             :    * = 2 z^(m-1) sum_{k >= 1} zeta(2k+2) * Z^(k+1) / (2k+2)..(2k+1+m)), where
    1906             :    * Z = -(z/2Pi)^2. Stop at 2k = (li - (m-1)*Lz - m) / dz, Lz = log2 |z| */
    1907             :   /* We cut the sum in two: small values of k first */
    1908          21 :   Z = gsqr(z); av = avma;
    1909          21 :   ksmall = get_k(dz, prec2nbits(prec));
    1910          21 :   constbern(ksmall);
    1911         469 :   for(k = 1; k < ksmall; k++)
    1912             :   {
    1913         469 :     GEN t = q = gdivgunextu(gmul(q,Z), 2*k+m); /* z^(2k+m+1)/(2k+m+1)! */
    1914         469 :     if (real) t = real_i(t);
    1915         469 :     t = gmul(t, gdivgu(bernfrac(2*k+2), 2*k+2)); /* - t * zeta(1-(2k+2)) */
    1916         469 :     s = gsub(s, t);
    1917         469 :     if (gexpo(t)  < li) return s;
    1918             :     /* large values ? */
    1919         448 :     if ((k & 0x1ff) == 0) gerepileall(av, 2, &s, &q);
    1920             :   }
    1921           0 :   if (DEBUGLEVEL>2) timer_printf(&T, "polylog: small k <= %ld", k);
    1922           0 :   Z = gneg(gsqr(gdiv(z, Pi2n(1,prec))));
    1923           0 :   q = gmul(gpowgs(z, m-1), gpowgs(Z, k+1)); /* Z^(k+1) * z^(m-1) */
    1924           0 :   S = gen_0; av = avma;
    1925           0 :   for(;; k++)
    1926           0 :   {
    1927           0 :     GEN t = q;
    1928             :     long b;
    1929           0 :     if (real) t = real_i(t);
    1930           0 :     b = prec + gexpo(t) / BITS_IN_LONG; /* decrease accuracy */
    1931           0 :     if (b == 2) break;
    1932             :     /* t * zeta(2k+2) / (2k+2)..(2k+1+m) */
    1933           0 :     t = gdiv(t, mulri(inv_szeta_euler(2*k+2, b),
    1934           0 :                       mulu_interval(2*k+2, 2*k+1+m)));
    1935           0 :     S = gadd(S, t); if (gexpo(t)  < li) break;
    1936           0 :     q = gmul(q, Z);
    1937           0 :     if ((k & 0x1ff) == 0) gerepileall(av, 2, &S, &q);
    1938             :   }
    1939           0 :   if (DEBUGLEVEL>2) timer_printf(&T, "polylog: large k <= %ld", k);
    1940           0 :   return gadd(s, gmul2n(S,1));
    1941             : }
    1942             : 
    1943             : static GEN
    1944          42 : Li1(GEN x, long prec) { return gneg(glog(gsubsg(1, x), prec)); }
    1945             : 
    1946             : static GEN
    1947         203 : polylog(long m, GEN x, long prec)
    1948             : {
    1949             :   long l, e, i, G, sx;
    1950             :   pari_sp av, av1;
    1951             :   GEN X, Xn, z, p1, p2, y, res;
    1952             : 
    1953         203 :   if (m < 0) pari_err_DOMAIN("polylog", "index", "<", gen_0, stoi(m));
    1954         203 :   if (!m) return mkfrac(gen_m1,gen_2);
    1955         203 :   if (gequal0(x)) return gcopy(x);
    1956         203 :   if (m==1) { av = avma; return gerepileupto(av, Li1(x, prec)); }
    1957             : 
    1958         168 :   l = precision(x);
    1959         168 :   if (!l) l = prec; else prec = l;
    1960         168 :   res = cgetc(l); av = avma;
    1961         168 :   x = gtofp(x, l+EXTRAPRECWORD);
    1962         168 :   e = gexpo(gnorm(x));
    1963         168 :   if (!e || e == -1) {
    1964          21 :     y = cxpolylog(m,x,prec);
    1965          21 :     set_avma(av); return affc_fixlg(y, res);
    1966             :   }
    1967         147 :   X = (e > 0)? ginv(x): x;
    1968         147 :   G = -prec2nbits(l);
    1969         147 :   av1 = avma;
    1970         147 :   y = Xn = X;
    1971         147 :   for (i=2; ; i++)
    1972             :   {
    1973       68159 :     Xn = gmul(X,Xn); p2 = gdiv(Xn,powuu(i,m));
    1974       68159 :     y = gadd(y,p2);
    1975       68159 :     if (gexpo(p2) <= G) break;
    1976             : 
    1977       68012 :     if (gc_needed(av1,1))
    1978             :     {
    1979           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"polylog");
    1980           0 :       gerepileall(av1,2, &y, &Xn);
    1981             :     }
    1982             :   }
    1983         147 :   if (e < 0) { set_avma(av); return affc_fixlg(y, res); }
    1984             : 
    1985          28 :   sx = gsigne(imag_i(x));
    1986          28 :   if (!sx)
    1987             :   {
    1988          28 :     if (m&1) sx = gsigne(gsub(gen_1, real_i(x)));
    1989          21 :     else     sx = - gsigne(real_i(x));
    1990             :   }
    1991          28 :   z = divri(mppi(l), mpfact(m-1)); setsigne(z, sx);
    1992          28 :   z = mkcomplex(gen_0, z);
    1993             : 
    1994          28 :   if (m == 2)
    1995             :   { /* same formula as below, written more efficiently */
    1996          21 :     y = gneg_i(y);
    1997          21 :     if (typ(x) == t_REAL && signe(x) < 0)
    1998           7 :       p1 = logr_abs(x);
    1999             :     else
    2000          14 :       p1 = gsub(glog(x,l), z);
    2001          21 :     p1 = gmul2n(gsqr(p1), -1); /* = (log(-x))^2 / 2 */
    2002             : 
    2003          21 :     p1 = gadd(p1, divru(sqrr(mppi(l)), 6));
    2004          21 :     p1 = gneg_i(p1);
    2005             :   }
    2006             :   else
    2007             :   {
    2008           7 :     GEN logx = glog(x,l), logx2 = gsqr(logx), vz = constzeta(m, l);
    2009           7 :     p1 = mkfrac(gen_m1,gen_2);
    2010          14 :     for (i = m-2; i >= 0; i -= 2)
    2011           7 :       p1 = gadd(gel(vz, m-i), gmul(p1, gdivgunextu(logx2, i+1)));
    2012           7 :     if (m&1) p1 = gmul(logx,p1); else y = gneg_i(y);
    2013           7 :     p1 = gadd(gmul2n(p1,1), gmul(z,gpowgs(logx,m-1)));
    2014           7 :     if (typ(x) == t_REAL && signe(x) < 0) p1 = real_i(p1);
    2015             :   }
    2016          28 :   y = gadd(y,p1);
    2017          28 :   set_avma(av); return affc_fixlg(y, res);
    2018             : }
    2019             : static GEN
    2020         119 : RIpolylog(long m, GEN x, long real, long prec)
    2021             : {
    2022         119 :   GEN y = polylog(m, x, prec);
    2023         119 :   return real? real_i(y): imag_i(y);
    2024             : }
    2025             : GEN
    2026          21 : dilog(GEN x, long prec) { return gpolylog(2, x, prec); }
    2027             : 
    2028             : /* x a floating point number, t_REAL or t_COMPLEX of t_REAL */
    2029             : static GEN
    2030          42 : logabs(GEN x)
    2031             : {
    2032             :   GEN y;
    2033          42 :   if (typ(x) == t_COMPLEX)
    2034             :   {
    2035           7 :     y = logr_abs( cxnorm(x) );
    2036           7 :     shiftr_inplace(y, -1);
    2037             :   } else
    2038          35 :     y = logr_abs(x);
    2039          42 :   return y;
    2040             : }
    2041             : 
    2042             : static GEN
    2043          21 : polylogD(long m, GEN x, long flag, long prec)
    2044             : {
    2045          21 :   long fl = 0, k, l, m2;
    2046             :   pari_sp av;
    2047             :   GEN p1, p2, y;
    2048             : 
    2049          21 :   if (gequal0(x)) return gcopy(x);
    2050          21 :   m2 = m&1;
    2051          21 :   if (gequal1(x) && m>=2) return m2? szeta(m,prec): gen_0;
    2052          21 :   av = avma; l = precision(x);
    2053          21 :   if (!l) { l = prec; x = gtofp(x,l); }
    2054          21 :   p1 = logabs(x);
    2055          21 :   if (signe(p1) > 0) { x = ginv(x); fl = !m2; } else setabssign(p1);
    2056             :   /* |x| <= 1, p1 = - log|x| >= 0 */
    2057          21 :   p2 = gen_1;
    2058          21 :   y = RIpolylog(m, x, m2, l);
    2059          84 :   for (k = 1; k < m; k++)
    2060             :   {
    2061          63 :     GEN t = RIpolylog(m-k, x, m2, l);
    2062          63 :     p2 = gdivgu(gmul(p2,p1), k); /* (-log|x|)^k / k! */
    2063          63 :     y = gadd(y, gmul(p2, t));
    2064             :   }
    2065          21 :   if (m2)
    2066             :   {
    2067          14 :     p1 = flag? gdivgs(p1, -2*m): gdivgs(logabs(gsubsg(1,x)), m);
    2068          14 :     y = gadd(y, gmul(p2, p1));
    2069             :   }
    2070          21 :   if (fl) y = gneg(y);
    2071          21 :   return gerepileupto(av, y);
    2072             : }
    2073             : 
    2074             : static GEN
    2075          14 : polylogP(long m, GEN x, long prec)
    2076             : {
    2077          14 :   long fl = 0, k, l, m2;
    2078             :   pari_sp av;
    2079             :   GEN p1,y;
    2080             : 
    2081          14 :   if (gequal0(x)) return gcopy(x);
    2082          14 :   m2 = m&1;
    2083          14 :   if (gequal1(x) && m>=2) return m2? szeta(m,prec): gen_0;
    2084          14 :   av = avma; l = precision(x);
    2085          14 :   if (!l) { l = prec; x = gtofp(x,l); }
    2086          14 :   p1 = logabs(x);
    2087          14 :   if (signe(p1) > 0) { x = ginv(x); fl = !m2; setsigne(p1, -1); }
    2088             :   /* |x| <= 1 */
    2089          14 :   y = RIpolylog(m, x, m2, l);
    2090          14 :   if (m==1)
    2091             :   {
    2092           7 :     shiftr_inplace(p1, -1); /* log |x| / 2 */
    2093           7 :     y = gadd(y, p1);
    2094             :   }
    2095             :   else
    2096             :   { /* m >= 2, \sum_{0 <= k <= m} 2^k B_k/k! (log |x|)^k Li_{m-k}(x),
    2097             :        with Li_0(x) := -1/2 */
    2098           7 :     GEN u, t = RIpolylog(m-1, x, m2, l);
    2099           7 :     u = gneg_i(p1); /* u = 2 B1 log |x| */
    2100           7 :     y = gadd(y, gmul(u, t));
    2101           7 :     if (m > 2)
    2102             :     {
    2103             :       GEN p2;
    2104           7 :       shiftr_inplace(p1, 1); /* 2log|x| <= 0 */
    2105           7 :       constbern(m>>1);
    2106           7 :       p1 = sqrr(p1);
    2107           7 :       p2 = shiftr(p1,-1);
    2108          21 :       for (k = 2; k < m; k += 2)
    2109             :       {
    2110          14 :         if (k > 2) p2 = gdivgunextu(gmul(p2,p1),k-1); /* 2^k/k! log^k |x|*/
    2111          14 :         t = RIpolylog(m-k, x, m2, l);
    2112          14 :         u = gmul(p2, bernfrac(k));
    2113          14 :         y = gadd(y, gmul(u, t));
    2114             :       }
    2115             :     }
    2116             :   }
    2117          14 :   if (fl) y = gneg(y);
    2118          14 :   return gerepileupto(av, y);
    2119             : }
    2120             : 
    2121             : static GEN
    2122         175 : gpolylog_i(void *E, GEN x, long prec)
    2123             : {
    2124         175 :   pari_sp av = avma;
    2125         175 :   long i, n, v, m = (long)E;
    2126             :   GEN a, y;
    2127             : 
    2128         175 :   if (m <= 0)
    2129             :   {
    2130          28 :     a = gmul(x, poleval(eulerianpol(-m, 0), x));
    2131          28 :     return gerepileupto(av, gdiv(a, gpowgs(gsubsg(1, x), 1-m)));
    2132             :   }
    2133         147 :   switch(typ(x))
    2134             :   {
    2135          84 :     case t_REAL: case t_COMPLEX: return polylog(m,x,prec);
    2136           7 :     case t_INTMOD: case t_PADIC: pari_err_IMPL( "padic polylogarithm");
    2137          56 :     default:
    2138          56 :       av = avma; if (!(y = toser_i(x))) break;
    2139          21 :       if (!m) { set_avma(av); return mkfrac(gen_m1,gen_2); }
    2140          21 :       if (m==1) return gerepileupto(av, Li1(y, prec));
    2141          21 :       if (gequal0(y)) return gerepilecopy(av, y);
    2142          21 :       v = valp(y);
    2143          21 :       if (v < 0) pari_err_DOMAIN("polylog","valuation", "<", gen_0, x);
    2144          14 :       if (v > 0) {
    2145           7 :         n = (lg(y)-3 + v) / v;
    2146           7 :         a = zeroser(varn(y), lg(y)-2);
    2147          35 :         for (i=n; i>=1; i--)
    2148          28 :           a = gmul(y, gadd(a, powis(utoipos(i),-m)));
    2149             :       } else { /* v == 0 */
    2150           7 :         long vy = varn(y);
    2151           7 :         GEN a0 = polcoef_i(y, 0, -1), t = gdiv(derivser(y), y);
    2152           7 :         a = Li1(y, prec);
    2153          14 :         for (i=2; i<=m; i++)
    2154           7 :           a = gadd(gpolylog(i, a0, prec), integ(gmul(t, a), vy));
    2155             :       }
    2156          14 :       return gerepileupto(av, a);
    2157             :   }
    2158          35 :   return trans_evalgen("polylog", E, gpolylog_i, x, prec);
    2159             : }
    2160             : GEN
    2161         133 : gpolylog(long m, GEN x, long prec) { return gpolylog_i((void*)m, x, prec); }
    2162             : 
    2163             : GEN
    2164         147 : polylog0(long m, GEN x, long flag, long prec)
    2165             : {
    2166         147 :   switch(flag)
    2167             :   {
    2168         105 :     case 0: return gpolylog(m,x,prec);
    2169          14 :     case 1: return polylogD(m,x,0,prec);
    2170           7 :     case 2: return polylogD(m,x,1,prec);
    2171          14 :     case 3: return polylogP(m,x,prec);
    2172           7 :     default: pari_err_FLAG("polylog");
    2173             :   }
    2174             :   return NULL; /* LCOV_EXCL_LINE */
    2175             : }
    2176             : 
    2177             : GEN
    2178       22411 : upper_to_cx(GEN x, long *prec)
    2179             : {
    2180       22411 :   long tx = typ(x), l;
    2181       22411 :   if (tx == t_QUAD) { x = quadtofp(x, *prec); tx = typ(x); }
    2182       22411 :   switch(tx)
    2183             :   {
    2184       22390 :     case t_COMPLEX:
    2185       22390 :       if (gsigne(gel(x,2)) > 0) break; /*fall through*/
    2186             :     case t_REAL: case t_INT: case t_FRAC:
    2187          14 :       pari_err_DOMAIN("modular function", "Im(argument)", "<=", gen_0, x);
    2188           7 :     default:
    2189           7 :       pari_err_TYPE("modular function", x);
    2190             :   }
    2191       22390 :   l = precision(x); if (l) *prec = l;
    2192       22390 :   return x;
    2193             : }
    2194             : 
    2195             : /* sqrt(3)/2 */
    2196             : static GEN
    2197        2189 : sqrt32(long prec) { GEN z = sqrtr_abs(utor(3,prec)); setexpo(z, -1); return z; }
    2198             : /* exp(i x), x = k pi/12 */
    2199             : static GEN
    2200        3532 : e12(ulong k, long prec)
    2201             : {
    2202             :   int s, sPi, sPiov2;
    2203             :   GEN z, t;
    2204        3532 :   k %= 24;
    2205        3532 :   if (!k) return gen_1;
    2206        3532 :   if (k == 12) return gen_m1;
    2207        3532 :   if (k >12) { s = 1; k = 24 - k; } else s = 0; /* x -> 2pi - x */
    2208        3532 :   if (k > 6) { sPi = 1; k = 12 - k; } else sPi = 0; /* x -> pi  - x */
    2209        3532 :   if (k > 3) { sPiov2 = 1; k = 6 - k; } else sPiov2 = 0; /* x -> pi/2 - x */
    2210        3532 :   z = cgetg(3, t_COMPLEX);
    2211        3532 :   switch(k)
    2212             :   {
    2213         630 :     case 0: gel(z,1) = icopy(gen_1); gel(z,2) = gen_0; break;
    2214         782 :     case 1: t = gmul2n(addrs(sqrt32(prec), 1), -1);
    2215         782 :       gel(z,1) = sqrtr(t);
    2216         782 :       gel(z,2) = gmul2n(invr(gel(z,1)), -2); break;
    2217             : 
    2218        1407 :     case 2: gel(z,1) = sqrt32(prec);
    2219        1407 :             gel(z,2) = real2n(-1, prec); break;
    2220             : 
    2221         713 :     case 3: gel(z,1) = sqrtr_abs(real2n(-1,prec));
    2222         713 :             gel(z,2) = rcopy(gel(z,1)); break;
    2223             :   }
    2224        3532 :   if (sPiov2) swap(gel(z,1), gel(z,2));
    2225        3532 :   if (sPi) togglesign(gel(z,1));
    2226        3532 :   if (s)   togglesign(gel(z,2));
    2227        3532 :   return z;
    2228             : }
    2229             : /* z a t_FRAC */
    2230             : static GEN
    2231       12548 : expIPifrac(GEN z, long prec)
    2232             : {
    2233       12548 :   GEN n = gel(z,1), d = gel(z,2);
    2234       12548 :   ulong r, q = uabsdivui_rem(12, d, &r);
    2235       12548 :   if (!r) return e12(q * umodiu(n, 24), prec); /* 12 | d */
    2236        9016 :   n = centermodii(n, shifti(d,1), d);
    2237        9016 :   return expIr(divri(mulri(mppi(prec), n), d));
    2238             : }
    2239             : /* exp(i Pi z), z a t_INT or t_FRAC */
    2240             : static GEN
    2241       10997 : expIPiQ(GEN z, long prec)
    2242             : {
    2243       10997 :   if (typ(z) == t_INT) return mpodd(z)? gen_m1: gen_1;
    2244        1958 :   return expIPifrac(z, prec);
    2245             : }
    2246             : 
    2247             : /* convert power of 2 t_REAL to rational */
    2248             : static GEN
    2249        1719 : real2nQ(GEN x)
    2250             : {
    2251        1719 :   long e = expo(x);
    2252             :   GEN z;
    2253        1719 :   if (e < 0)
    2254         426 :     z = mkfrac(signe(x) < 0? gen_m1: gen_1, int2n(-e));
    2255             :   else
    2256             :   {
    2257        1293 :     z = int2n(e);
    2258        1293 :     if (signe(x) < 0) togglesign_safe(&z);
    2259             :   }
    2260        1719 :   return z;
    2261             : }
    2262             : /* x a real number */
    2263             : GEN
    2264      180775 : expIPiR(GEN x, long prec)
    2265             : {
    2266      180775 :   if (typ(x) == t_REAL && absrnz_equal2n(x)) x = real2nQ(x);
    2267      180775 :   switch(typ(x))
    2268             :   {
    2269        3036 :     case t_INT:  return mpodd(x)? gen_m1: gen_1;
    2270        1941 :     case t_FRAC: return expIPifrac(x, prec);
    2271             :   }
    2272      175798 :   return expIr(mulrr(mppi(prec), x));
    2273             : }
    2274             : /* z a t_COMPLEX */
    2275             : GEN
    2276      232740 : expIPiC(GEN z, long prec)
    2277             : {
    2278             :   GEN pi, r, x, y;
    2279      232740 :   if (typ(z) != t_COMPLEX) return expIPiR(z, prec);
    2280       52252 :   x = gel(z,1);
    2281       52252 :   y = gel(z,2); if (gequal0(y)) return expIPiR(x, prec);
    2282       51965 :   pi = mppi(prec);
    2283       51965 :   r = gmul(pi, y); togglesign(r); r = mpexp(r); /* exp(-pi y) */
    2284       51965 :   if (typ(x) == t_REAL && absrnz_equal2n(x)) x = real2nQ(x);
    2285       51965 :   switch(typ(x))
    2286             :   {
    2287        4550 :     case t_INT: if (mpodd(x)) togglesign(r);
    2288        4550 :                 return r;
    2289        8649 :     case t_FRAC: return gmul(r, expIPifrac(x, prec));
    2290             :   }
    2291       38766 :   return gmul(r, expIr(mulrr(pi, x)));
    2292             : }
    2293             : /* exp(I x y), more efficient for x in R, y pure imaginary */
    2294             : GEN
    2295         147 : expIxy(GEN x, GEN y, long prec) { return gexp(gmul(x, mulcxI(y)), prec); }
    2296             : 
    2297             : static GEN
    2298       13549 : qq(GEN x, long prec)
    2299             : {
    2300       13549 :   long tx = typ(x);
    2301             :   GEN y;
    2302             : 
    2303       13549 :   if (is_scalar_t(tx))
    2304             :   {
    2305       13507 :     if (tx == t_PADIC) return x;
    2306       13493 :     x = upper_to_cx(x, &prec);
    2307       13479 :     return cxtoreal(expIPiC(gmul2n(x,1), prec)); /* e(x) */
    2308             :   }
    2309          42 :   if (! ( y = toser_i(x)) ) pari_err_TYPE("modular function", x);
    2310          42 :   return y;
    2311             : }
    2312             : 
    2313             : /* return (y * X^d) + x. Assume d > 0, x != 0, valp(x) = 0 */
    2314             : static GEN
    2315          21 : ser_addmulXn(GEN y, GEN x, long d)
    2316             : {
    2317          21 :   long i, lx, ly, l = valp(y) + d; /* > 0 */
    2318             :   GEN z;
    2319             : 
    2320          21 :   lx = lg(x);
    2321          21 :   ly = lg(y) + l; if (lx < ly) ly = lx;
    2322          21 :   if (l > lx-2) return gcopy(x);
    2323          21 :   z = cgetg(ly,t_SER);
    2324          77 :   for (i=2; i<=l+1; i++) gel(z,i) = gel(x,i);
    2325          70 :   for (   ; i < ly; i++) gel(z,i) = gadd(gel(x,i),gel(y,i-l));
    2326          21 :   z[1] = x[1]; return z;
    2327             : }
    2328             : 
    2329             : /* q a t_POL s.t. q(0) != 0, v > 0, Q = x^v*q; return \prod_i (1-Q^i) */
    2330             : static GEN
    2331          28 : RgXn_eta(GEN q, long v, long lim)
    2332             : {
    2333          28 :   pari_sp av = avma;
    2334             :   GEN qn, ps, y;
    2335             :   ulong vps, vqn, n;
    2336             : 
    2337          28 :   if (!degpol(q) && isint1(gel(q,2))) return eta_ZXn(v, lim+v);
    2338           7 :   y = qn = ps = pol_1(0);
    2339           7 :   vps = vqn = 0;
    2340           7 :   for(n = 0;; n++)
    2341           7 :   { /* qn = q^n,  ps = (-1)^n q^(n(3n+1)/2),
    2342             :      * vps, vqn valuation of ps, qn HERE */
    2343          14 :     pari_sp av2 = avma;
    2344          14 :     ulong vt = vps + 2*vqn + v; /* valuation of t at END of loop body */
    2345             :     long k1, k2;
    2346             :     GEN t;
    2347          14 :     vqn += v; vps = vt + vqn; /* valuation of qn, ps at END of body */
    2348          14 :     k1 = lim + v - vt + 1;
    2349          14 :     k2 = k1 - vqn; /* = lim + v - vps + 1 */
    2350          14 :     if (k1 <= 0) break;
    2351          14 :     t = RgX_mul(q, RgX_sqr(qn));
    2352          14 :     t = RgXn_red_shallow(t, k1);
    2353          14 :     t = RgX_mul(ps,t);
    2354          14 :     t = RgXn_red_shallow(t, k1);
    2355          14 :     t = RgX_neg(t); /* t = (-1)^(n+1) q^(n(3n+1)/2 + 2n+1) */
    2356          14 :     t = gerepileupto(av2, t);
    2357          14 :     y = RgX_addmulXn_shallow(t, y, vt);
    2358          14 :     if (k2 <= 0) break;
    2359             : 
    2360           7 :     qn = RgX_mul(qn,q);
    2361           7 :     ps = RgX_mul(t,qn);
    2362           7 :     ps = RgXn_red_shallow(ps, k2);
    2363           7 :     y = RgX_addmulXn_shallow(ps, y, vps);
    2364             : 
    2365           7 :     if (gc_needed(av,1))
    2366             :     {
    2367           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"eta, n = %ld", n);
    2368           0 :       gerepileall(av, 3, &y, &qn, &ps);
    2369             :     }
    2370             :   }
    2371           7 :   return y;
    2372             : }
    2373             : 
    2374             : static GEN
    2375       16455 : inteta(GEN q)
    2376             : {
    2377       16455 :   long tx = typ(q);
    2378             :   GEN ps, qn, y;
    2379             : 
    2380       16455 :   y = gen_1; qn = gen_1; ps = gen_1;
    2381       16455 :   if (tx==t_PADIC)
    2382             :   {
    2383          28 :     if (valp(q) <= 0) pari_err_DOMAIN("eta", "v_p(q)", "<=",gen_0,q);
    2384             :     for(;;)
    2385          56 :     {
    2386          77 :       GEN t = gneg_i(gmul(ps,gmul(q,gsqr(qn))));
    2387          77 :       y = gadd(y,t); qn = gmul(qn,q); ps = gmul(t,qn);
    2388          77 :       t = y;
    2389          77 :       y = gadd(y,ps); if (gequal(t,y)) return y;
    2390             :     }
    2391             :   }
    2392             : 
    2393       16427 :   if (tx == t_SER)
    2394             :   {
    2395             :     ulong vps, vqn;
    2396          42 :     long l = lg(q), v, n;
    2397             :     pari_sp av;
    2398             : 
    2399          42 :     v = valp(q); /* handle valuation separately to avoid overflow */
    2400          42 :     if (v <= 0) pari_err_DOMAIN("eta", "v_p(q)", "<=",gen_0,q);
    2401          35 :     y = ser2pol_i(q, l); /* t_SER inefficient when input has low degree */
    2402          35 :     n = degpol(y);
    2403          35 :     if (n <= (l>>2))
    2404             :     {
    2405          28 :       GEN z = RgXn_eta(y, v, l-2);
    2406          28 :       setvarn(z, varn(y)); return RgX_to_ser(z, l+v);
    2407             :     }
    2408           7 :     q = leafcopy(q); av = avma;
    2409           7 :     setvalp(q, 0);
    2410           7 :     y = scalarser(gen_1, varn(q), l+v);
    2411           7 :     vps = vqn = 0;
    2412           7 :     for(n = 0;; n++)
    2413           7 :     { /* qn = q^n,  ps = (-1)^n q^(n(3n+1)/2) */
    2414          14 :       ulong vt = vps + 2*vqn + v;
    2415             :       long k;
    2416             :       GEN t;
    2417          14 :       t = gneg_i(gmul(ps,gmul(q,gsqr(qn))));
    2418             :       /* t = (-1)^(n+1) q^(n(3n+1)/2 + 2n+1) */
    2419          14 :       y = ser_addmulXn(t, y, vt);
    2420          14 :       vqn += v; vps = vt + vqn;
    2421          14 :       k = l+v - vps; if (k <= 2) return y;
    2422             : 
    2423           7 :       qn = gmul(qn,q); ps = gmul(t,qn);
    2424           7 :       y = ser_addmulXn(ps, y, vps);
    2425           7 :       setlg(q, k);
    2426           7 :       setlg(qn, k);
    2427           7 :       setlg(ps, k);
    2428           7 :       if (gc_needed(av,3))
    2429             :       {
    2430           0 :         if(DEBUGMEM>1) pari_warn(warnmem,"eta");
    2431           0 :         gerepileall(av, 3, &y, &qn, &ps);
    2432             :       }
    2433             :     }
    2434             :   }
    2435             :   {
    2436       16385 :     long l = -prec2nbits(precision(q));
    2437       16385 :     pari_sp av = avma;
    2438             : 
    2439             :     for(;;)
    2440       53573 :     {
    2441       69958 :       GEN t = gneg_i(gmul(ps,gmul(q,gsqr(qn))));
    2442             :       /* qn = q^n
    2443             :        * ps = (-1)^n q^(n(3n+1)/2)
    2444             :        * t = (-1)^(n+1) q^(n(3n+1)/2 + 2n+1) */
    2445       69958 :       y = gadd(y,t); qn = gmul(qn,q); ps = gmul(t,qn);
    2446       69958 :       y = gadd(y,ps);
    2447       69958 :       if (gexpo(ps)-gexpo(y) < l) return y;
    2448       53573 :       if (gc_needed(av,3))
    2449             :       {
    2450           0 :         if(DEBUGMEM>1) pari_warn(warnmem,"eta");
    2451           0 :         gerepileall(av, 3, &y, &qn, &ps);
    2452             :       }
    2453             :     }
    2454             :   }
    2455             : }
    2456             : 
    2457             : GEN
    2458          77 : eta(GEN x, long prec)
    2459             : {
    2460          77 :   pari_sp av = avma;
    2461          77 :   GEN z = inteta( qq(x,prec) );
    2462          49 :   if (typ(z) == t_SER) return gerepilecopy(av, z);
    2463          14 :   return gerepileupto(av, z);
    2464             : }
    2465             : 
    2466             : /* s(h,k) = sum(n = 1, k-1, (n/k)*(frac(h*n/k) - 1/2))
    2467             :  * Knuth's algorithm. h integer, k integer > 0, (h,k) = 1 */
    2468             : GEN
    2469        6258 : sumdedekind_coprime(GEN h, GEN k)
    2470             : {
    2471        6258 :   pari_sp av = avma;
    2472             :   GEN s2, s1, p, pp;
    2473             :   long s;
    2474        6258 :   if (lgefint(k) == 3 && uel(k,2) <= (2*(ulong)LONG_MAX) / 3)
    2475             :   {
    2476        6251 :     ulong kk = k[2], hh = umodiu(h, kk);
    2477             :     long s1, s2;
    2478             :     GEN v;
    2479        6251 :     if (signe(k) < 0) { k = negi(k); hh = Fl_neg(hh, kk); }
    2480        6251 :     v = u_sumdedekind_coprime(hh, kk);
    2481        6251 :     s1 = v[1]; s2 = v[2];
    2482        6251 :     return gerepileupto(av, gdiv(addis(mulis(k,s1), s2), muluu(12, kk)));
    2483             :   }
    2484           7 :   s = 1;
    2485           7 :   s1 = gen_0; p = gen_1; pp = gen_0;
    2486           7 :   s2 = h = modii(h, k);
    2487          35 :   while (signe(h)) {
    2488          28 :     GEN r, nexth, a = dvmdii(k, h, &nexth);
    2489          28 :     if (is_pm1(h)) s2 = s == 1? addii(s2, p): subii(s2, p);
    2490          28 :     s1 = s == 1? addii(s1, a): subii(s1, a);
    2491          28 :     s = -s;
    2492          28 :     k = h; h = nexth;
    2493          28 :     r = addii(mulii(a,p), pp); pp = p; p = r;
    2494             :   }
    2495             :   /* at this point p = original k */
    2496           7 :   if (s == -1) s1 = subiu(s1, 3);
    2497           7 :   return gerepileupto(av, gdiv(addii(mulii(p,s1), s2), muliu(p,12)));
    2498             : }
    2499             : /* as above, for ulong arguments.
    2500             :  * k integer > 0, 0 <= h < k, (h,k) = 1. Returns [s1,s2] such that
    2501             :  * s(h,k) = (s2 + k s1) / (12k). Requires max(h + k/2, k) < LONG_MAX
    2502             :  * to avoid overflow, in particular k <= LONG_MAX * 2/3 is fine */
    2503             : GEN
    2504        6251 : u_sumdedekind_coprime(long h, long k)
    2505             : {
    2506        6251 :   long s = 1, s1 = 0, s2 = h, p = 1, pp = 0;
    2507       11424 :   while (h) {
    2508        5173 :     long r, nexth = k % h, a = k / h; /* a >= 1, a >= 2 if h = 1 */
    2509        5173 :     if (h == 1) s2 += p * s; /* occurs exactly once, last step */
    2510        5173 :     s1 += a * s;
    2511        5173 :     s = -s;
    2512        5173 :     k = h; h = nexth;
    2513        5173 :     r = a*p + pp; pp = p; p = r; /* p >= pp >= 0 */
    2514             :   }
    2515             :   /* in the above computation, p increases from 1 to original k,
    2516             :    * -k/2 <= s2 <= h + k/2, and |s1| <= k */
    2517        6251 :   if (s < 0) s1 -= 3; /* |s1| <= k+3 ? */
    2518             :   /* But in fact, |s2 + p s1| <= k^2 + 1/2 - 3k; if (s < 0), we have
    2519             :    * |s2| <= k/2 and it follows that |s1| < k here as well */
    2520             :   /* p = k; s(h,k) = (s2 + p s1)/12p. */
    2521        6251 :   return mkvecsmall2(s1, s2);
    2522             : }
    2523             : GEN
    2524          28 : sumdedekind(GEN h, GEN k)
    2525             : {
    2526          28 :   pari_sp av = avma;
    2527             :   GEN d;
    2528          28 :   if (typ(h) != t_INT) pari_err_TYPE("sumdedekind",h);
    2529          28 :   if (typ(k) != t_INT) pari_err_TYPE("sumdedekind",k);
    2530          28 :   d = gcdii(h,k);
    2531          28 :   if (!is_pm1(d))
    2532             :   {
    2533           7 :     h = diviiexact(h, d);
    2534           7 :     k = diviiexact(k, d);
    2535             :   }
    2536          28 :   return gerepileupto(av, sumdedekind_coprime(h,k));
    2537             : }
    2538             : 
    2539             : /* eta(x); assume Im x >> 0 (e.g. x in SL2's standard fundamental domain) */
    2540             : static GEN
    2541       17297 : eta_reduced(GEN x, long prec)
    2542             : {
    2543       17297 :   GEN z = expIPiC(gdivgu(x, 12), prec); /* e(x/24) */
    2544       17297 :   if (24 * gexpo(z) >= -prec2nbits(prec))
    2545       16364 :     z = gmul(z, inteta( gpowgs(z,24) ));
    2546       17297 :   return z;
    2547             : }
    2548             : 
    2549             : /* x = U.z (flag = 1), or x = U^(-1).z (flag = 0)
    2550             :  * Return [s,t] such that eta(z) = eta(x) * sqrt(s) * exp(I Pi t) */
    2551             : static GEN
    2552       17311 : eta_correction(GEN x, GEN U, long flag)
    2553             : {
    2554             :   GEN a,b,c,d, s,t;
    2555             :   long sc;
    2556       17311 :   a = gcoeff(U,1,1);
    2557       17311 :   b = gcoeff(U,1,2);
    2558       17311 :   c = gcoeff(U,2,1);
    2559       17311 :   d = gcoeff(U,2,2);
    2560             :   /* replace U by U^(-1) */
    2561       17311 :   if (flag) {
    2562        8939 :     swap(a,d);
    2563        8939 :     togglesign_safe(&b);
    2564        8939 :     togglesign_safe(&c);
    2565             :   }
    2566       17311 :   sc = signe(c);
    2567       17311 :   if (!sc) {
    2568       11081 :     if (signe(d) < 0) togglesign_safe(&b);
    2569       11081 :     s = gen_1;
    2570       11081 :     t = uutoQ(umodiu(b, 24), 12);
    2571             :   } else {
    2572        6230 :     if (sc < 0) {
    2573        1722 :       togglesign_safe(&a);
    2574        1722 :       togglesign_safe(&b);
    2575        1722 :       togglesign_safe(&c);
    2576        1722 :       togglesign_safe(&d);
    2577             :     } /* now c > 0 */
    2578        6230 :     s = mulcxmI(gadd(gmul(c,x), d));
    2579        6230 :     t = gadd(gdiv(addii(a,d),muliu(c,12)), sumdedekind_coprime(negi(d),c));
    2580             :     /* correction : exp(I Pi (((a+d)/12c) + s(-d,c)) ) sqrt(-i(cx+d))  */
    2581             :   }
    2582       17311 :   return mkvec2(s, t);
    2583             : }
    2584             : 
    2585             : /* returns the true value of eta(x) for Im(x) > 0, using reduction to
    2586             :  * standard fundamental domain */
    2587             : GEN
    2588        8869 : trueeta(GEN x, long prec)
    2589             : {
    2590        8869 :   pari_sp av = avma;
    2591             :   GEN U, st, s, t;
    2592             : 
    2593        8869 :   if (!is_scalar_t(typ(x))) pari_err_TYPE("trueeta",x);
    2594        8869 :   x = upper_to_cx(x, &prec);
    2595        8869 :   x = cxredsl2(x, &U);
    2596        8869 :   st = eta_correction(x, U, 1);
    2597        8869 :   x = eta_reduced(x, prec);
    2598        8869 :   s = gel(st, 1);
    2599        8869 :   t = gel(st, 2);
    2600        8869 :   x = gmul(x, expIPiQ(t, prec));
    2601        8869 :   if (s != gen_1) x = gmul(x, gsqrt(s, prec));
    2602        8869 :   return gerepileupto(av, x);
    2603             : }
    2604             : 
    2605             : GEN
    2606          77 : eta0(GEN x, long flag,long prec)
    2607          77 : { return flag? trueeta(x,prec): eta(x,prec); }
    2608             : 
    2609             : /* eta(q) = 1 + \sum_{n>0} (-1)^n * (q^(n(3n-1)/2) + q^(n(3n+1)/2)) */
    2610             : static GEN
    2611           7 : ser_eta(long prec)
    2612             : {
    2613           7 :   GEN e = cgetg(prec+2, t_SER), ed = e+2;
    2614             :   long n, j;
    2615           7 :   e[1] = evalsigne(1)|_evalvalp(0)|evalvarn(0);
    2616           7 :   gel(ed,0) = gen_1;
    2617         483 :   for (n = 1; n < prec; n++) gel(ed,n) = gen_0;
    2618          49 :   for (n = 1, j = 0; n < prec; n++)
    2619             :   {
    2620             :     GEN s;
    2621          49 :     j += 3*n-2; /* = n*(3*n-1) / 2 */;
    2622          49 :     if (j >= prec) break;
    2623          42 :     s = odd(n)? gen_m1: gen_1;
    2624          42 :     gel(ed, j) = s;
    2625          42 :     if (j+n >= prec) break;
    2626          42 :     gel(ed, j+n) = s;
    2627             :   }
    2628           7 :   return e;
    2629             : }
    2630             : 
    2631             : static GEN
    2632         476 : coeffEu(GEN fa)
    2633             : {
    2634         476 :   pari_sp av = avma;
    2635         476 :   return gerepileuptoint(av, mului(65520, usumdivk_fact(fa,11)));
    2636             : }
    2637             : /* E12 = 1 + q*E/691 */
    2638             : static GEN
    2639           7 : ser_E(long prec)
    2640             : {
    2641           7 :   GEN e = cgetg(prec+2, t_SER), ed = e+2;
    2642           7 :   GEN F = vecfactoru_i(2, prec); /* F[n] = factoru(n+1) */
    2643             :   long n;
    2644           7 :   e[1] = evalsigne(1)|_evalvalp(0)|evalvarn(0);
    2645           7 :   gel(ed,0) = utoipos(65520);
    2646         483 :   for (n = 1; n < prec; n++) gel(ed,n) = coeffEu(gel(F,n));
    2647           7 :   return e;
    2648             : }
    2649             : /* j = E12/Delta + 432000/691, E12 = 1 + q*E/691 */
    2650             : static GEN
    2651           7 : ser_j2(long prec, long v)
    2652             : {
    2653           7 :   pari_sp av = avma;
    2654           7 :   GEN iD = gpowgs(ginv(ser_eta(prec)), 24); /* q/Delta */
    2655           7 :   GEN J = gmul(ser_E(prec), iD);
    2656           7 :   setvalp(iD,-1); /* now 1/Delta */
    2657           7 :   J = gadd(gdivgu(J, 691), iD);
    2658           7 :   J = gerepileupto(av, J);
    2659           7 :   if (prec > 1) gel(J,3) = utoipos(744);
    2660           7 :   setvarn(J,v); return J;
    2661             : }
    2662             : 
    2663             : /* j(q) = \sum_{n >= -1} c(n)q^n,
    2664             :  * \sum_{n = -1}^{N-1} c(n) (-10n \sigma_3(N-n) + 21 \sigma_5(N-n))
    2665             :  * = c(N) (N+1)/24 */
    2666             : static GEN
    2667          14 : ser_j(long prec, long v)
    2668             : {
    2669             :   GEN j, J, S3, S5, F;
    2670             :   long i, n;
    2671          14 :   if (prec > 64) return ser_j2(prec, v);
    2672           7 :   S3 = cgetg(prec+1, t_VEC);
    2673           7 :   S5 = cgetg(prec+1,t_VEC);
    2674           7 :   F = vecfactoru_i(1, prec);
    2675          35 :   for (n = 1; n <= prec; n++)
    2676             :   {
    2677          28 :     GEN fa = gel(F,n);
    2678          28 :     gel(S3,n) = mului(10, usumdivk_fact(fa,3));
    2679          28 :     gel(S5,n) = mului(21, usumdivk_fact(fa,5));
    2680             :   }
    2681           7 :   J = cgetg(prec+2, t_SER),
    2682           7 :   J[1] = evalvarn(v)|evalsigne(1)|evalvalp(-1);
    2683           7 :   j = J+3;
    2684           7 :   gel(j,-1) = gen_1;
    2685           7 :   gel(j,0) = utoipos(744);
    2686           7 :   gel(j,1) = utoipos(196884);
    2687          21 :   for (n = 2; n < prec; n++)
    2688             :   {
    2689          14 :     pari_sp av = avma;
    2690          14 :     GEN c, s3 = gel(S3,n+1), s5 = gel(S5,n+1);
    2691          14 :     c = addii(s3, s5);
    2692          49 :     for (i = 0; i < n; i++)
    2693             :     {
    2694          35 :       s3 = gel(S3,n-i); s5 = gel(S5,n-i);
    2695          35 :       c = addii(c, mulii(gel(j,i), subii(s5, mului(i,s3))));
    2696             :     }
    2697          14 :     gel(j,n) = gerepileuptoint(av, diviuexact(muliu(c,24), n+1));
    2698             :   }
    2699           7 :   return J;
    2700             : }
    2701             : 
    2702             : GEN
    2703          42 : jell(GEN x, long prec)
    2704             : {
    2705          42 :   long tx = typ(x);
    2706          42 :   pari_sp av = avma;
    2707             :   GEN q, h, U;
    2708             : 
    2709          42 :   if (!is_scalar_t(tx))
    2710             :   {
    2711             :     long v;
    2712          21 :     if (gequalX(x)) return ser_j(precdl, varn(x));
    2713          21 :     q = toser_i(x); if (!q) pari_err_TYPE("ellj",x);
    2714          14 :     v = fetch_var_higher();
    2715          14 :     h = ser_j(lg(q)-2, v);
    2716          14 :     h = gsubst(h, v, q);
    2717          14 :     delete_var(); return gerepileupto(av, h);
    2718             :   }
    2719          21 :   if (tx == t_PADIC)
    2720             :   {
    2721           7 :     GEN p2, p1 = gdiv(inteta(gsqr(x)), inteta(x));
    2722           7 :     p1 = gmul2n(gsqr(p1),1);
    2723           7 :     p1 = gmul(x,gpowgs(p1,12));
    2724           7 :     p2 = gaddsg(768,gadd(gsqr(p1),gdivsg(4096,p1)));
    2725           7 :     p1 = gmulsg(48,p1);
    2726           7 :     return gerepileupto(av, gadd(p2,p1));
    2727             :   }
    2728             :   /* Let h = Delta(2x) / Delta(x), then j(x) = (1 + 256h)^3 / h */
    2729          14 :   x = upper_to_cx(x, &prec);
    2730           7 :   x = cxredsl2(x, &U); /* forget about Ua : j has weight 0 */
    2731             :   { /* cf eta_reduced, raised to power 24
    2732             :      * Compute
    2733             :      *   t = (inteta(q(2x)) / inteta(q(x))) ^ 24;
    2734             :      * then
    2735             :      *   h = t * (q(2x) / q(x) = t * q(x);
    2736             :      * but inteta(q) costly and useless if expo(q) << 1  => inteta(q) = 1.
    2737             :      * log_2 ( exp(-2Pi Im tau) ) < -prec2nbits(prec)
    2738             :      * <=> Im tau > prec2nbits(prec) * log(2) / 2Pi */
    2739           7 :     long C = (long)prec2nbits_mul(prec, M_LN2/(2*M_PI));
    2740           7 :     q = expIPiC(gmul2n(x,1), prec); /* e(x) */
    2741           7 :     if (gcmpgs(gel(x,2), C) > 0) /* eta(q(x)) = 1 : no need to compute q(2x) */
    2742           0 :       h = q;
    2743             :     else
    2744             :     {
    2745           7 :       GEN t = gdiv(inteta(gsqr(q)), inteta(q));
    2746           7 :       h = gmul(q, gpowgs(t, 24));
    2747             :     }
    2748             :   }
    2749             :   /* real_1 important ! gaddgs(, 1) could increase the accuracy ! */
    2750           7 :   return gerepileupto(av, gdiv(gpowgs(gadd(gmul2n(h,8), real_1(prec)), 3), h));
    2751             : }
    2752             : 
    2753             : static GEN
    2754        8372 : to_form(GEN a, GEN w, GEN C)
    2755        8372 : { return mkvec3(a, w, diviiexact(C, a)); }
    2756             : static GEN
    2757        8372 : form_to_quad(GEN f, GEN sqrtD)
    2758             : {
    2759        8372 :   long a = itos(gel(f,1)), a2 = a << 1;
    2760        8372 :   GEN b = gel(f,2);
    2761        8372 :   return mkcomplex(gdivgs(b, -a2), gdivgs(sqrtD, a2));
    2762             : }
    2763             : static GEN
    2764        8372 : eta_form(GEN f, GEN sqrtD, GEN *s_t, long prec)
    2765             : {
    2766        8372 :   GEN U, t = form_to_quad(redimagsl2(f, &U), sqrtD);
    2767        8372 :   *s_t = eta_correction(t, U, 0);
    2768        8372 :   return eta_reduced(t, prec);
    2769             : }
    2770             : 
    2771             : /* eta(t/p)eta(t/q) / (eta(t)eta(t/pq)), t = (-w + sqrt(D)) / 2a */
    2772             : GEN
    2773        2093 : double_eta_quotient(GEN a, GEN w, GEN D, long p, long q, GEN pq, GEN sqrtD)
    2774             : {
    2775        2093 :   GEN C = shifti(subii(sqri(w), D), -2);
    2776             :   GEN d, t, z, zp, zq, zpq, s_t, s_tp, s_tpq, s, sp, spq;
    2777        2093 :   long prec = realprec(sqrtD);
    2778             : 
    2779        2093 :   z = eta_form(to_form(a, w, C), sqrtD, &s_t, prec);
    2780        2093 :   s = gel(s_t, 1);
    2781        2093 :   zp = eta_form(to_form(mului(p, a), w, C), sqrtD, &s_tp, prec);
    2782        2093 :   sp = gel(s_tp, 1);
    2783        2093 :   zpq = eta_form(to_form(mulii(pq, a), w, C), sqrtD, &s_tpq, prec);
    2784        2093 :   spq = gel(s_tpq, 1);
    2785        2093 :   if (p == q) {
    2786           0 :     z = gdiv(gsqr(zp), gmul(z, zpq));
    2787           0 :     t = gsub(gmul2n(gel(s_tp,2), 1),
    2788           0 :              gadd(gel(s_t,2), gel(s_tpq,2)));
    2789           0 :     if (sp != gen_1) z = gmul(z, sp);
    2790             :   } else {
    2791             :     GEN s_tq, sq;
    2792        2093 :     zq = eta_form(to_form(mului(q, a), w, C), sqrtD, &s_tq, prec);
    2793        2093 :     sq = gel(s_tq, 1);
    2794        2093 :     z = gdiv(gmul(zp, zq), gmul(z, zpq));
    2795        2093 :     t = gsub(gadd(gel(s_tp,2), gel(s_tq,2)),
    2796        2093 :              gadd(gel(s_t,2), gel(s_tpq,2)));
    2797        2093 :     if (sp != gen_1) z = gmul(z, gsqrt(sp, prec));
    2798        2093 :     if (sq != gen_1) z = gmul(z, gsqrt(sq, prec));
    2799             :   }
    2800        2093 :   d = NULL;
    2801        2093 :   if (s != gen_1) d = gsqrt(s, prec);
    2802        2093 :   if (spq != gen_1) {
    2803        2065 :     GEN x = gsqrt(spq, prec);
    2804        2065 :     d = d? gmul(d, x): x;
    2805             :   }
    2806        2093 :   if (d) z = gdiv(z, d);
    2807        2093 :   return gmul(z, expIPiQ(t, prec));
    2808             : }
    2809             : 
    2810             : typedef struct { GEN u; long v, t; } cxanalyze_t;
    2811             : 
    2812             : /* Check whether a t_COMPLEX, t_REAL or t_INT z != 0 can be written as
    2813             :  * z = u * 2^(v/2) * exp(I Pi/4 t), u > 0, v = 0,1 and -3 <= t <= 4.
    2814             :  * Allow z t_INT/t_REAL to simplify handling of eta_correction() output */
    2815             : static int
    2816          70 : cxanalyze(cxanalyze_t *T, GEN z)
    2817             : {
    2818             :   GEN a, b;
    2819             :   long ta, tb;
    2820             : 
    2821          70 :   T->v = 0;
    2822          70 :   if (is_intreal_t(typ(z)))
    2823             :   {
    2824          63 :     T->u = mpabs_shallow(z);
    2825          63 :     T->t = signe(z) < 0? 4: 0;
    2826          63 :     return 1;
    2827             :   }
    2828           7 :   a = gel(z,1); ta = typ(a);
    2829           7 :   b = gel(z,2); tb = typ(b);
    2830             : 
    2831           7 :   T->t = 0;
    2832           7 :   if (ta == t_INT && !signe(a))
    2833             :   {
    2834           0 :     T->u = R_abs_shallow(b);
    2835           0 :     T->t = gsigne(b) < 0? -2: 2;
    2836           0 :     return 1;
    2837             :   }
    2838           7 :   if (tb == t_INT && !signe(b))
    2839             :   {
    2840           0 :     T->u = R_abs_shallow(a);
    2841           0 :     T->t = gsigne(a) < 0? 4: 0;
    2842           0 :     return 1;
    2843             :   }
    2844           7 :   if (ta != tb || ta == t_REAL) { T->u = z; return 0; }
    2845             :   /* a,b both non zero, both t_INT or t_FRAC */
    2846           7 :   if (ta == t_INT)
    2847             :   {
    2848           7 :     if (!absequalii(a, b)) return 0;
    2849           7 :     T->u = absi_shallow(a);
    2850           7 :     T->v = 1;
    2851           7 :     if (signe(a) == signe(b))
    2852           0 :     { T->t = signe(a) < 0? -3: 1; }
    2853             :     else
    2854           7 :     { T->t = signe(a) < 0? 3: -1; }
    2855             :   }
    2856             :   else
    2857             :   {
    2858           0 :     if (!absequalii(gel(a,2), gel(b,2)) || !absequalii(gel(a,1),gel(b,1)))
    2859           0 :       return 0;
    2860           0 :     T->u = absfrac_shallow(a);
    2861           0 :     T->v = 1;
    2862           0 :     a = gel(a,1);
    2863           0 :     b = gel(b,1);
    2864           0 :     if (signe(a) == signe(b))
    2865           0 :     { T->t = signe(a) < 0? -3: 1; }
    2866             :     else
    2867           0 :     { T->t = signe(a) < 0? 3: -1; }
    2868             :   }
    2869           7 :   return 1;
    2870             : }
    2871             : 
    2872             : /* z * sqrt(st_b) / sqrt(st_a) exp(I Pi (t + t0)). Assume that
    2873             :  * sqrt2 = gsqrt(gen_2, prec) or NULL */
    2874             : static GEN
    2875          35 : apply_eta_correction(GEN z, GEN st_a, GEN st_b, GEN t0, GEN sqrt2, long prec)
    2876             : {
    2877          35 :   GEN t, s_a = gel(st_a, 1), s_b = gel(st_b, 1);
    2878             :   cxanalyze_t Ta, Tb;
    2879             :   int ca, cb;
    2880             : 
    2881          35 :   t = gsub(gel(st_b,2), gel(st_a,2));
    2882          35 :   if (t0 != gen_0) t = gadd(t, t0);
    2883          35 :   ca = cxanalyze(&Ta, s_a);
    2884          35 :   cb = cxanalyze(&Tb, s_b);
    2885          35 :   if (ca || cb)
    2886          35 :   { /* compute sqrt(s_b) / sqrt(s_a) in a more efficient way:
    2887             :      * sb = ub sqrt(2)^vb exp(i Pi/4 tb) */
    2888          35 :     GEN u = gdiv(Tb.u,Ta.u);
    2889          35 :     switch(Tb.v - Ta.v)
    2890             :     {
    2891           0 :       case -1: u = gmul2n(u,-1); /* fall through: write 1/sqrt2 = sqrt2/2 */
    2892           7 :       case 1: u = gmul(u, sqrt2? sqrt2: sqrtr_abs(real2n(1, prec)));
    2893             :     }
    2894          35 :     if (!isint1(u)) z = gmul(z, gsqrt(u, prec));
    2895          35 :     t = gadd(t, gmul2n(stoi(Tb.t - Ta.t), -3));
    2896             :   }
    2897             :   else
    2898             :   {
    2899           0 :     z = gmul(z, gsqrt(s_b, prec));
    2900           0 :     z = gdiv(z, gsqrt(s_a, prec));
    2901             :   }
    2902          35 :   return gmul(z, expIPiQ(t, prec));
    2903             : }
    2904             : 
    2905             : /* sqrt(2) eta(2x) / eta(x) */
    2906             : GEN
    2907           7 : weberf2(GEN x, long prec)
    2908             : {
    2909           7 :   pari_sp av = avma;
    2910             :   GEN z, sqrt2, a,b, Ua,Ub, st_a,st_b;
    2911             : 
    2912           7 :   x = upper_to_cx(x, &prec);
    2913           7 :   a = cxredsl2(x, &Ua);
    2914           7 :   b = cxredsl2(gmul2n(x,1), &Ub);
    2915           7 :   if (gequal(a,b)) /* not infrequent */
    2916           0 :     z = gen_1;
    2917             :   else
    2918           7 :     z = gdiv(eta_reduced(b,prec), eta_reduced(a,prec));
    2919           7 :   st_a = eta_correction(a, Ua, 1);
    2920           7 :   st_b = eta_correction(b, Ub, 1);
    2921           7 :   sqrt2 = sqrtr_abs(real2n(1, prec));
    2922           7 :   z = apply_eta_correction(z, st_a, st_b, gen_0, sqrt2, prec);
    2923           7 :   return gerepileupto(av, gmul(z, sqrt2));
    2924             : }
    2925             : 
    2926             : /* eta(x/2) / eta(x) */
    2927             : GEN
    2928          14 : weberf1(GEN x, long prec)
    2929             : {
    2930          14 :   pari_sp av = avma;
    2931             :   GEN z, a,b, Ua,Ub, st_a,st_b;
    2932             : 
    2933          14 :   x = upper_to_cx(x, &prec);
    2934          14 :   a = cxredsl2(x, &Ua);
    2935          14 :   b = cxredsl2(gmul2n(x,-1), &Ub);
    2936          14 :   if (gequal(a,b)) /* not infrequent */
    2937           0 :     z = gen_1;
    2938             :   else
    2939          14 :     z = gdiv(eta_reduced(b,prec), eta_reduced(a,prec));
    2940          14 :   st_a = eta_correction(a, Ua, 1);
    2941          14 :   st_b = eta_correction(b, Ub, 1);
    2942          14 :   z = apply_eta_correction(z, st_a, st_b, gen_0, NULL, prec);
    2943          14 :   return gerepileupto(av, z);
    2944             : }
    2945             : /* exp(-I*Pi/24) * eta((x+1)/2) / eta(x) */
    2946             : GEN
    2947          14 : weberf(GEN x, long prec)
    2948             : {
    2949          14 :   pari_sp av = avma;
    2950             :   GEN z, t0, a,b, Ua,Ub, st_a,st_b;
    2951          14 :   x = upper_to_cx(x, &prec);
    2952          14 :   a = cxredsl2(x, &Ua);
    2953          14 :   b = cxredsl2(gmul2n(gaddgs(x,1),-1), &Ub);
    2954          14 :   if (gequal(a,b)) /* not infrequent */
    2955           7 :     z = gen_1;
    2956             :   else
    2957           7 :     z = gdiv(eta_reduced(b,prec), eta_reduced(a,prec));
    2958          14 :   st_a = eta_correction(a, Ua, 1);
    2959          14 :   st_b = eta_correction(b, Ub, 1);
    2960          14 :   t0 = mkfrac(gen_m1, utoipos(24));
    2961          14 :   z = apply_eta_correction(z, st_a, st_b, t0, NULL, prec);
    2962          14 :   if (typ(z) == t_COMPLEX && isexactzero(real_i(x)))
    2963           0 :     z = gerepilecopy(av, gel(z,1));
    2964             :   else
    2965          14 :     z = gerepileupto(av, z);
    2966          14 :   return z;
    2967             : }
    2968             : GEN
    2969          35 : weber0(GEN x, long flag,long prec)
    2970             : {
    2971          35 :   switch(flag)
    2972             :   {
    2973          14 :     case 0: return weberf(x,prec);
    2974          14 :     case 1: return weberf1(x,prec);
    2975           7 :     case 2: return weberf2(x,prec);
    2976           0 :     default: pari_err_FLAG("weber");
    2977             :   }
    2978             :   return NULL; /* LCOV_EXCL_LINE */
    2979             : }
    2980             : 
    2981             : /* check |q| < 1 */
    2982             : static GEN
    2983          21 : check_unit_disc(const char *fun, GEN q, long prec)
    2984             : {
    2985          21 :   GEN Q = gtofp(q, prec), Qlow;
    2986          21 :   Qlow = (prec > LOWDEFAULTPREC)? gtofp(Q,LOWDEFAULTPREC): Q;
    2987          21 :   if (gcmp(gnorm(Qlow), gen_1) >= 0)
    2988           0 :     pari_err_DOMAIN(fun, "abs(q)", ">=", gen_1, q);
    2989          21 :   return Q;
    2990             : }
    2991             : 
    2992             : GEN
    2993          14 : theta(GEN q, GEN z, long prec)
    2994             : {
    2995             :   long l, n;
    2996          14 :   pari_sp av = avma, av2;
    2997             :   GEN s, c, snz, cnz, s2z, c2z, ps, qn, y, zy, ps2, k, zold;
    2998             : 
    2999          14 :   l = precision(q);
    3000          14 :   n = precision(z); if (n && n < l) l = n;
    3001          14 :   if (l) prec = l;
    3002          14 :   z = gtofp(z, prec);
    3003          14 :   q = check_unit_disc("theta", q, prec);
    3004          14 :   zold = NULL; /* gcc -Wall */
    3005          14 :   zy = imag_i(z);
    3006          14 :   if (gequal0(zy)) k = gen_0;
    3007             :   else
    3008             :   {
    3009           7 :     GEN lq = glog(q,prec); k = roundr(divrr(zy, real_i(lq)));
    3010           7 :     if (signe(k)) { zold = z; z = gadd(z, mulcxmI(gmul(lq,k))); }
    3011             :   }
    3012          14 :   qn = gen_1;
    3013          14 :   ps2 = gsqr(q);
    3014          14 :   ps = gneg_i(ps2);
    3015          14 :   gsincos(z, &s, &c, prec);
    3016          14 :   s2z = gmul2n(gmul(s,c),1); /* sin 2z */
    3017          14 :   c2z = gsubgs(gmul2n(gsqr(c),1), 1); /* cos 2z */
    3018          14 :   snz = s;
    3019          14 :   cnz = c; y = s;
    3020          14 :   av2 = avma;
    3021          14 :   for (n = 3;; n += 2)
    3022         259 :   {
    3023             :     long e;
    3024         273 :     s = gadd(gmul(snz, c2z), gmul(cnz,s2z));
    3025         273 :     qn = gmul(qn,ps);
    3026         273 :     y = gadd(y, gmul(qn, s));
    3027         273 :     e = gexpo(s); if (e < 0) e = 0;
    3028         273 :     if (gexpo(qn) + e < -prec2nbits(prec)) break;
    3029             : 
    3030         259 :     ps = gmul(ps,ps2);
    3031         259 :     c = gsub(gmul(cnz, c2z), gmul(snz,s2z));
    3032         259 :     snz = s; /* sin nz */
    3033         259 :     cnz = c; /* cos nz */
    3034         259 :     if (gc_needed(av,2))
    3035             :     {
    3036           0 :       if (DEBUGMEM>1) pari_warn(warnmem,"theta (n = %ld)", n);
    3037           0 :       gerepileall(av2, 5, &snz, &cnz, &ps, &qn, &y);
    3038             :     }
    3039             :   }
    3040          14 :   if (signe(k))
    3041             :   {
    3042           7 :     y = gmul(y, gmul(powgi(q,sqri(k)),
    3043             :                      gexp(gmul(mulcxI(zold),shifti(k,1)), prec)));
    3044           7 :     if (mod2(k)) y = gneg_i(y);
    3045             :   }
    3046          14 :   return gerepileupto(av, gmul(y, gmul2n(gsqrt(gsqrt(q,prec),prec),1)));
    3047             : }
    3048             : 
    3049             : GEN
    3050           7 : thetanullk(GEN q, long k, long prec)
    3051             : {
    3052             :   long l, n;
    3053           7 :   pari_sp av = avma;
    3054             :   GEN p1, ps, qn, y, ps2;
    3055             : 
    3056           7 :   if (k < 0)
    3057           0 :     pari_err_DOMAIN("thetanullk", "k", "<", gen_0, stoi(k));
    3058           7 :   l = precision(q);
    3059           7 :   if (l) prec = l;
    3060           7 :   q = check_unit_disc("thetanullk", q, prec);
    3061             : 
    3062           7 :   if (!odd(k)) { set_avma(av); return gen_0; }
    3063           7 :   qn = gen_1;
    3064           7 :   ps2 = gsqr(q);
    3065           7 :   ps = gneg_i(ps2);
    3066           7 :   y = gen_1;
    3067           7 :   for (n = 3;; n += 2)
    3068         280 :   {
    3069             :     GEN t;
    3070         287 :     qn = gmul(qn,ps);
    3071         287 :     ps = gmul(ps,ps2);
    3072         287 :     t = gmul(qn, powuu(n, k)); y = gadd(y, t);
    3073         287 :     if (gexpo(t) < -prec2nbits(prec)) break;
    3074             :   }
    3075           7 :   p1 = gmul2n(gsqrt(gsqrt(q,prec),prec),1);
    3076           7 :   if (k&2) y = gneg_i(y);
    3077           7 :   return gerepileupto(av, gmul(p1, y));
    3078             : }
    3079             : 
    3080             : /* q2 = q^2 */
    3081             : static GEN
    3082       13472 : vecthetanullk_loop(GEN q2, long k, long prec)
    3083             : {
    3084       13472 :   GEN ps, qn = gen_1, y = const_vec(k, gen_1);
    3085       13472 :   pari_sp av = avma;
    3086       13472 :   const long bit = prec2nbits(prec);
    3087             :   long i, n;
    3088             : 
    3089       13472 :   if (gexpo(q2) < -2*bit) return y;
    3090       13472 :   ps = gneg_i(q2);
    3091       13472 :   for (n = 3;; n += 2)
    3092       95595 :   {
    3093      109067 :     GEN t = NULL/*-Wall*/, P = utoipos(n), N2 = sqru(n);
    3094      109067 :     qn = gmul(qn,ps);
    3095      109067 :     ps = gmul(ps,q2);
    3096      327201 :     for (i = 1; i <= k; i++)
    3097             :     {
    3098      218134 :       t = gmul(qn, P); gel(y,i) = gadd(gel(y,i), t);
    3099      218134 :       P = mulii(P, N2);
    3100             :     }
    3101      109067 :     if (gexpo(t) < -bit) return y;
    3102       95595 :     if (gc_needed(av,2))
    3103             :     {
    3104           0 :       if (DEBUGMEM>1) pari_warn(warnmem,"vecthetanullk_loop, n = %ld",n);
    3105           0 :       gerepileall(av, 3, &qn, &ps, &y);
    3106             :     }
    3107             :   }
    3108             : }
    3109             : /* [d^i theta/dz^i(q, 0), i = 1, 3, .., 2*k - 1] */
    3110             : GEN
    3111           0 : vecthetanullk(GEN q, long k, long prec)
    3112             : {
    3113           0 :   long i, l = precision(q);
    3114           0 :   pari_sp av = avma;
    3115             :   GEN p1, y;
    3116             : 
    3117           0 :   if (l) prec = l;
    3118           0 :   q = check_unit_disc("vecthetanullk", q, prec);
    3119           0 :   y = vecthetanullk_loop(gsqr(q), k, prec);
    3120           0 :   p1 = gmul2n(gsqrt(gsqrt(q,prec),prec),1);
    3121           0 :   for (i = 2; i <= k; i += 2) gel(y,i) = gneg_i(gel(y,i));
    3122           0 :   return gerepileupto(av, gmul(p1, y));
    3123             : }
    3124             : 
    3125             : /* [d^i theta/dz^i(q, 0), i = 1, 3, .., 2*k - 1], q = exp(2iPi tau) */
    3126             : GEN
    3127           0 : vecthetanullk_tau(GEN tau, long k, long prec)
    3128             : {
    3129           0 :   long i, l = precision(tau);
    3130           0 :   pari_sp av = avma;
    3131             :   GEN q4, y;
    3132             : 
    3133           0 :   if (l) prec = l;
    3134           0 :   if (typ(tau) != t_COMPLEX || gsigne(gel(tau,2)) <= 0)
    3135           0 :     pari_err_DOMAIN("vecthetanullk_tau", "imag(tau)", "<=", gen_0, tau);
    3136           0 :   q4 = expIPiC(gmul2n(tau,-1), prec); /* q^(1/4) */
    3137           0 :   y = vecthetanullk_loop(gpowgs(q4,8), k, prec);
    3138           0 :   for (i = 2; i <= k; i += 2) gel(y,i) = gneg_i(gel(y,i));
    3139           0 :   return gerepileupto(av, gmul(gmul2n(q4,1), y));
    3140             : }
    3141             : 
    3142             : /* Return E_k(tau). Slow if tau is not in standard fundamental domain */
    3143             : GEN
    3144       13741 : cxEk(GEN tau, long k, long prec)
    3145             : {
    3146             :   pari_sp av;
    3147             :   GEN q, y, qn;
    3148       13741 :   long n, b, l = precision(tau);
    3149             : 
    3150       13741 :   if (l) prec = l;
    3151       13741 :   b = bit_accuracy(prec);
    3152             :   /* sum n^(k-1) x^n <= x(1 + (k!-1)x) / (1-x)^k (cf Eulerian polynomials)
    3153             :    * S = \sum_{n > 0} n^(k-1) |q^n/(1-q^n)| <= x(1+(k!-1)x) / (1-x)^(k+1),
    3154             :    * where x = |q| = exp(-2Pi Im(tau)) < 1. Neglegt 2/zeta(1-k) * S if
    3155             :    * (2Pi)^k/(k-1)! x < 2^(-b-1) and k! x < 1. Use log2((2Pi)^k/(k-1)!) < 10 */
    3156       13741 :   if (gcmpgs(imag_i(tau), (M_LN2 / (2*M_PI)) * (b+1+10)) > 0)
    3157          17 :     return real_1(prec);
    3158       13724 :   if (k == 2)
    3159             :   { /* -theta^(3)(tau/2) / theta^(1)(tau/2). Assume that Im tau > 0 */
    3160       13472 :     y = vecthetanullk_loop(qq(tau,prec), 2, prec);
    3161       13472 :     return gdiv(gel(y,2), gel(y,1));
    3162             :   }
    3163         252 :   q = cxtoreal(expIPiC(gneg(gmul2n(tau,1)), prec));
    3164         252 :   av = avma; y = gen_0; qn = q;
    3165         252 :   for(n = 1;; n++)
    3166        2605 :   { /* compute y := sum_{n>0} n^(k-1) / (q^n-1) */
    3167        2857 :     GEN t = gdiv(powuu(n,k-1), gsubgs(qn,1));
    3168        2857 :     if (gequal0(t) || gexpo(t) <= - prec2nbits(prec) - 5) break;
    3169        2605 :     y = gadd(y, t);
    3170        2605 :     qn = gmul(q,qn);
    3171        2605 :     if (gc_needed(av,2))
    3172             :     {
    3173           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"elleisnum");
    3174           0 :       gerepileall(av, 2, &y,&qn);
    3175             :     }
    3176             :   }
    3177         252 :   return gadd(gen_1, gmul(y, gdiv(gen_2, szeta(1-k, prec))));
    3178             : }

Generated by: LCOV version 1.13