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.12.1 lcov report (development 25406-bf255ab81b) Lines: 1903 2035 93.5 %
Date: 2020-06-04 05:59:24 Functions: 131 134 97.8 %
Legend: Lines: hit not hit

          Line data    Source code
       1             : /* Copyright (C) 2000  The PARI group.
       2             : 
       3             : This file is part of the PARI/GP package.
       4             : 
       5             : PARI/GP is free software; you can redistribute it and/or modify it under the
       6             : terms of the GNU General Public License as published by the Free Software
       7             : Foundation. It is distributed in the hope that it will be useful, but WITHOUT
       8             : ANY WARRANTY WHATSOEVER.
       9             : 
      10             : Check the License for details. You should have received a copy of it, along
      11             : with the package; see the file 'COPYING'. If not, write to the Free Software
      12             : Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. */
      13             : 
      14             : /********************************************************************/
      15             : /**                                                                **/
      16             : /**                   TRANSCENDENTAL FONCTIONS                     **/
      17             : /**                          (part 3)                              **/
      18             : /**                                                                **/
      19             : /********************************************************************/
      20             : #include "pari.h"
      21             : #include "paripriv.h"
      22             : 
      23             : #define HALF_E 1.3591409 /* exp(1) / 2 */
      24             : static const long EXTRAPREC = DEFAULTPREC-2;
      25             : 
      26             : /***********************************************************************/
      27             : /**                                                                   **/
      28             : /**                       BESSEL FUNCTIONS                            **/
      29             : /**                                                                   **/
      30             : /***********************************************************************/
      31             : 
      32             : static GEN
      33       11060 : _abs(GEN x)
      34       11060 : { return gabs(gtofp(x,LOWDEFAULTPREC), LOWDEFAULTPREC); }
      35             : /* can we use asymptotic expansion ? */
      36             : static int
      37        9135 : bessel_asymp(GEN z, long bit)
      38        9135 : { return gcmpgs(_abs(z), (bit+4)/2) >= 0; }
      39             : 
      40             : /* Region I: 0 < Arg z <= Pi, II: -Pi < Arg z <= 0 */
      41             : static int
      42         532 : regI(GEN z)
      43             : {
      44         532 :   long s = gsigne(imag_i(z));
      45         532 :   return (s > 0 || (s == 0 && gsigne(real_i(z)) < 0)) ? 1 : 2;
      46             : }
      47             : /* Region 1: Re(z) >= 0, 2: Re(z) < 0, Im(z) >= 0, 3: Re(z) < 0, Im(z) < 0 */
      48             : static int
      49         931 : regJ(GEN z)
      50             : {
      51         931 :   if (gsigne(real_i(z)) >= 0) return 1;
      52         336 :   return gsigne(imag_i(z)) >= 0 ? 2 : 3;
      53             : }
      54             : 
      55             : /* Hankel's expansions:
      56             :  * a_k(n) = \prod_{0 <= j < k} (4n^2 - (2j+1)^2)
      57             :  * C(k)[n,z] = a_k(n) / (k! (8 z)^k)
      58             :  * A(z)  = exp(-z) sum_{k >= 0} C(k)
      59             :  * A(-z) = exp(z) sum_{k >= 0} (-1)^k C(k)
      60             :  * J_n(z) ~ [1] (A(z/i) / r + A(-z/i) r) / sqrt(2Pi z)
      61             :  *          [2] (A(z/i) r^3 + A(-z/i) r) / sqrt(2Pi z)
      62             :  *          [3] (A(z/i) / r + A(-z/i) / r^3) / sqrt(2Pi z)
      63             :  * Y_n(z) ~ [1] i(A(z/i) / r + A(-z/i) r) / sqrt(2Pi z)
      64             :  *          [2] i(A(z/i) (r^3-2/r) + A(-z/i) r) / sqrt(2Pi z)
      65             :  *          [3] i(-A(z/i)/r + A(-z/i)(2r-1/r^3)) / sqrt(2Pi z)
      66             :  * K_n(z) ~ A(z) Pi / sqrt(2 Pi z)
      67             :  * I_n(z) ~ [I] (A(-z) + r^2 A(z)) / sqrt(2 Pi z)
      68             :  *          [II](A(-z) + r^(-2) A(z)) / sqrt(2 Pi z) */
      69             : 
      70             : /* set [A(z), A(-z), exp((2*nu+1)*I*Pi/4)] */
      71             : static void
      72        1925 : hankel_ABr(GEN *pA, GEN *pB, GEN *pr, GEN n, GEN z, long bit)
      73             : {
      74        1925 :   GEN E, P, C, Q = gen_0, zi = ginv(gmul2n(z, 3));
      75        1925 :   GEN K = gaddgs(_abs(n), 1), n2 = gmul2n(gsqr(n),2);
      76        1925 :   long prec = nbits2prec(bit), B = bit + 4, m;
      77             : 
      78        1925 :   P = C = real_1_bit(bit);
      79        1925 :   for (m = 1;; m += 2)
      80             :   {
      81       18466 :     C = gmul(C, gdivgs(gmul(gsub(n2, sqru(2*m - 1)), zi), m));
      82       18466 :     Q = gadd(Q, C);
      83       18466 :     C = gmul(C, gdivgs(gmul(gsub(n2, sqru(2*m + 1)), zi), m + 1));
      84       18466 :     P = gadd(P, C);
      85       18466 :     if (gexpo(C) < -B && gcmpgs(K, m) <= 0) break;
      86             :   }
      87        1925 :   E = gexp(z, prec);
      88        1925 :   *pA = gdiv(gadd(P, Q), E);
      89        1925 :   *pB = gmul(gsub(P, Q), E);
      90        1925 :   *pr = gexp(mulcxI(gmul(gaddgs(gmul2n(n,1), 1), Pi2n(-2, prec))), prec);
      91        1925 : }
      92             : 
      93             : /* sqrt(2*Pi*z) */
      94             : static GEN
      95        1925 : sqz(GEN z, long bit)
      96             : {
      97        1925 :   long prec = nbits2prec(bit);
      98        1925 :   return gsqrt(gmul(Pi2n(1, prec), z), prec);
      99             : }
     100             : 
     101             : static GEN
     102         462 : besskasymp(GEN nu, GEN z, long bit)
     103             : {
     104             :   GEN A, B, r;
     105         462 :   long prec = nbits2prec(bit);
     106         462 :   hankel_ABr(&A,&B,&r, nu, z, bit);
     107         462 :   return gdiv(gmul(A, mppi(prec)), sqz(z, bit));
     108             : }
     109             : 
     110             : static GEN
     111         532 : bessiasymp(GEN nu, GEN z, long bit)
     112             : {
     113             :   GEN A, B, r, R, r2;
     114         532 :   hankel_ABr(&A,&B,&r, nu, z, bit);
     115         532 :   r2 = gsqr(r);
     116         532 :   R = regI(z) == 1 ? gmul(A, r2) : gdiv(A, r2);
     117         532 :   return gdiv(gadd(B, R), sqz(z, bit));
     118             : }
     119             : 
     120             : static GEN
     121         469 : bessjasymp(GEN nu, GEN z, long bit)
     122             : {
     123             :   GEN A, B, r, R;
     124         469 :   long reg = regJ(z);
     125         469 :   hankel_ABr(&A,&B,&r, nu, mulcxmI(z), bit);
     126         469 :   if (reg == 1) R = gadd(gdiv(A, r), gmul(B, r));
     127         168 :   else if (reg == 2) R = gadd(gmul(A, gpowgs(r, 3)), gmul(B, r));
     128          56 :   else R = gadd(gdiv(A, r), gdiv(B, gpowgs(r, 3)));
     129         469 :   return gdiv(R, sqz(z, bit));
     130             : }
     131             : 
     132             : static GEN
     133         462 : bessyasymp(GEN nu, GEN z, long bit)
     134             : {
     135             :   GEN A, B, r, R;
     136         462 :   long reg = regJ(z);
     137         462 :   hankel_ABr(&A,&B,&r, nu, mulcxmI(z), bit);
     138         462 :   if (reg == 1) R = gsub(gmul(B, r), gdiv(A, r));
     139         168 :   else if (reg == 2)
     140         112 :     R = gadd(gmul(A, gsub(gpowgs(r, 3), gmul2n(ginv(r), 1))), gmul(B, r));
     141             :   else
     142          56 :     R = gsub(gmul(B, gsub(gmul2n(r, 1), ginv(gpowgs(r, 3)))), gdiv(A, r));
     143         462 :   return gdiv(mulcxI(R), sqz(z, bit));
     144             : }
     145             : 
     146             : /* n! sum_{k=0}^m Z^k / (k!*(k+n)!), with Z := (-1)^flag*z^2/4 */
     147             : static GEN
     148        4137 : _jbessel(GEN n, GEN z, long flag, long m)
     149             : {
     150             :   pari_sp av;
     151             :   GEN Z,s;
     152             :   long k;
     153             : 
     154        4137 :   Z = gmul2n(gsqr(z),-2); if (flag & 1) Z = gneg(Z);
     155        4137 :   if (typ(z) == t_SER)
     156             :   {
     157         210 :     long v = valp(z);
     158         210 :     if (v < 0) pari_err_DOMAIN("besselj","valuation", "<", gen_0, z);
     159         203 :     k = lg(Z)-2 - v;
     160         203 :     if (v == 0) pari_err_IMPL("besselj around a!=0");
     161         203 :     if (k <= 0) return scalarser(gen_1, varn(z), 2*v);
     162         203 :     setlg(Z, k+2);
     163             :   }
     164        4130 :   s = gen_1;
     165        4130 :   av = avma;
     166      187684 :   for (k=m; k>=1; k--)
     167             :   {
     168      183554 :     s = gaddsg(1, gdiv(gmul(Z,s), gmulgs(gaddgs(n,k),k)));
     169      183554 :     if (gc_needed(av,1))
     170             :     {
     171           0 :       if (DEBUGMEM>1) pari_warn(warnmem,"besselj");
     172           0 :       s = gerepileupto(av, s);
     173             :     }
     174             :   }
     175        4130 :   return s;
     176             : }
     177             : 
     178             : /* max(2, L * approximate solution to x log x = B) */
     179             : static long
     180        6629 : bessel_get_lim(double B, double L)
     181        6629 : { return maxss(2, L * exp(dbllambertW0(B))); }
     182             : 
     183             : static GEN jbesselintern(GEN n, GEN z, long flag, long prec);
     184             : static GEN kbesselintern(GEN n, GEN z, long flag, long prec);
     185             : static GEN
     186         140 : jbesselvec(GEN n, GEN x, long fl, long prec)
     187         308 : { pari_APPLY_same(jbesselintern(n,gel(x,i),fl,prec)) }
     188             : static GEN
     189          35 : jbesselhvec(GEN n, GEN x, long prec)
     190          77 : { pari_APPLY_same(jbesselh(n,gel(x,i),prec)) }
     191             : static GEN
     192          14 : polylogvec(long m, GEN x, long prec)
     193          42 : { pari_APPLY_same(gpolylog(m,gel(x,i),prec)) }
     194             : static GEN
     195         140 : kbesselvec(GEN n, GEN x, long fl, long prec)
     196         308 : { pari_APPLY_same(kbesselintern(n,gel(x,i),fl,prec)) }
     197             : 
     198             : /* flag = 0: I, flag = 1 : J */
     199             : static GEN
     200        5446 : jbesselintern(GEN n, GEN z, long flag, long prec)
     201             : {
     202             :   long i, ki;
     203        5446 :   pari_sp av = avma;
     204             :   GEN y;
     205             : 
     206        5446 :   switch(typ(z))
     207             :   {
     208           0 :     case t_QUAD: z = gtofp(z, prec); /* don't bother */
     209        5040 :     case t_INT: case t_FRAC: case t_REAL: case t_COMPLEX:
     210             :     {
     211        5040 :       int flz0 = gequal0(z);
     212             :       long lim, k, precnew, bit;
     213             :       GEN p1, p2;
     214             :       double az, L;
     215             : 
     216        5040 :       i = precision(z); if (i) prec = i;
     217        5040 :       if (flz0 && gequal0(n)) return real_1(prec);
     218        5040 :       bit = prec2nbits(prec);
     219        5040 :       if (bessel_asymp(z, bit))
     220             :       {
     221        1001 :         GEN R = (flag == 0)? bessiasymp(n, z, bit): bessjasymp(n, z, bit);
     222        1001 :         if (typ(R) == t_COMPLEX && isexactzero(imag_i(n))
     223         315 :                                 && gsigne(real_i(z)) > 0
     224         161 :                                 && isexactzero(imag_i(z))) R = gcopy(gel(R,1));
     225        1001 :         return gerepileupto(av, R);
     226             :       }
     227        4039 :       p2 = gpow(gmul2n(z,-1),n,prec);
     228        4011 :       p2 = gdiv(p2, ggamma(gaddgs(n,1),prec));
     229        4011 :       if (flz0) return gerepileupto(av, p2);
     230        4011 :       az = dblmodulus(z); L = HALF_E * az;
     231        4011 :       precnew = prec;
     232        4011 :       if (az >= 1.0) precnew += 1 + nbits2extraprec((long)(az/M_LN2));
     233        4011 :       if (issmall(n,&ki)) {
     234        3094 :         k = labs(ki);
     235        3094 :         n = utoi(k);
     236             :       } else {
     237         833 :         i = precision(n);
     238         833 :         if (i && i < precnew) n = gtofp(n,precnew);
     239             :       }
     240        3927 :       z = gtofp(z,precnew);
     241        3927 :       lim = bessel_get_lim(prec2nbits_mul(prec,M_LN2/2) / L, L);
     242        3927 :       p1 = gprec_wtrunc(_jbessel(n,z,flag,lim), prec);
     243        3927 :       return gerepileupto(av, gmul(p2,p1));
     244             :     }
     245             : 
     246         112 :     case t_VEC: case t_COL: case t_MAT:
     247         112 :       return jbesselvec(n, z, flag, prec);
     248             : 
     249          28 :     case t_POLMOD:
     250          28 :       y = jbesselvec(n, polmod_to_embed(z, prec), flag, prec);
     251          28 :       return gerepileupto(av,y);
     252             : 
     253          14 :     case t_PADIC: pari_err_IMPL("p-adic jbessel function");
     254         252 :     default:
     255         252 :       if (!(y = toser_i(z))) break;
     256         238 :       if (issmall(n,&ki)) n = utoi(labs(ki));
     257         210 :       return gerepileupto(av, _jbessel(n,y,flag,lg(y)-2));
     258             :   }
     259          14 :   pari_err_TYPE("jbessel",z);
     260             :   return NULL; /* LCOV_EXCL_LINE */
     261             : }
     262             : GEN
     263        1176 : jbessel(GEN n, GEN z, long prec) { return jbesselintern(n,z,1,prec); }
     264             : GEN
     265         847 : ibessel(GEN n, GEN z, long prec) { return jbesselintern(n,z,0,prec); }
     266             : 
     267             : /* k > 0 */
     268             : static GEN
     269         119 : _jbesselh(long k, GEN z, long prec)
     270             : {
     271         119 :   GEN s, c, p0, p1, zinv = ginv(z);
     272             :   long i;
     273             : 
     274         119 :   gsincos(z,&s,&c,prec);
     275         119 :   p1 = gmul(zinv,s);
     276         119 :   p0 = p1; p1 = gmul(zinv,gsub(p0,c));
     277        1134 :   for (i = 2; i <= k; i++)
     278             :   {
     279        1015 :     GEN p2 = gsub(gmul(gmulsg(2*i-1,zinv), p1), p0);
     280        1015 :     p0 = p1; p1 = p2;
     281             :   }
     282         119 :   return p1;
     283             : }
     284             : 
     285             : /* J_{n+1/2}(z) */
     286             : GEN
     287         315 : jbesselh(GEN n, GEN z, long prec)
     288             : {
     289             :   long k, i;
     290             :   pari_sp av;
     291             :   GEN y;
     292             : 
     293         315 :   if (typ(n)!=t_INT) pari_err_TYPE("jbesselh",n);
     294         203 :   k = itos(n);
     295         203 :   if (k < 0) return jbessel(gadd(ghalf,n), z, prec);
     296             : 
     297         203 :   switch(typ(z))
     298             :   {
     299           0 :     case t_QUAD: z = gtofp(z, prec); /* don't bother */
     300         133 :     case t_INT: case t_FRAC: case t_REAL: case t_COMPLEX:
     301             :     {
     302             :       long pr;
     303             :       GEN p1;
     304         133 :       if (gequal0(z))
     305             :       {
     306           7 :         av = avma;
     307           7 :         p1 = gmul(gsqrt(gdiv(z,mppi(prec)),prec),gpowgs(z,k));
     308           7 :         p1 = gdiv(p1, mulu_interval(k+1, 2*k+1)); /* x k! / (2k+1)! */
     309           7 :         return gerepileupto(av, gmul2n(p1,2*k));
     310             :       }
     311         126 :       if ( (pr = precision(z)) ) prec = pr;
     312         126 :       if (bessel_asymp(z, prec2nbits(prec)))
     313           7 :         return jbessel(gadd(ghalf,n), z, prec);
     314         119 :       y = cgetc(prec); av = avma;
     315         119 :       p1 = gsqrt(gdiv(z, Pi2n(-1,prec)), prec);
     316         119 :       if (!k)
     317          21 :         p1 = gmul(p1, gsinc(z, prec));
     318             :       else
     319             :       {
     320          98 :         long bits = BITS_IN_LONG + 2*k * (log2(k) -  dbllog2(z));
     321          98 :         if (bits > 0)
     322             :         {
     323          98 :           prec += nbits2extraprec(bits);
     324          98 :           if (pr) z = gtofp(z, prec);
     325             :         }
     326          98 :         p1 = gmul(p1, _jbesselh(k,z,prec));
     327             :       }
     328         119 :       set_avma(av); return affc_fixlg(p1, y);
     329             :     }
     330             : 
     331          28 :     case t_VEC: case t_COL: case t_MAT:
     332          28 :       return jbesselhvec(n, z, prec);
     333             : 
     334           7 :     case t_POLMOD:
     335           7 :       av = avma;
     336           7 :       return gerepileupto(av, jbesselhvec(n, polmod_to_embed(z,prec), prec));
     337             : 
     338           0 :     case t_PADIC: pari_err_IMPL("p-adic jbesselh function");
     339          35 :     default:
     340             :     {
     341             :       long t, v;
     342          35 :       av = avma; if (!(y = toser_i(z))) break;
     343          35 :       if (gequal0(y)) return gerepileupto(av, gpowgs(y,k));
     344          35 :       v = valp(y);
     345          35 :       if (v < 0) pari_err_DOMAIN("besseljh","valuation", "<", gen_0, y);
     346          28 :       t = lg(y)-2;
     347          28 :       if (v) y = sertoser(y, t + (2*k+1)*v);
     348          28 :       if (!k)
     349           7 :         y = gsinc(y,prec);
     350             :       else
     351             :       {
     352          21 :         GEN T, a = _jbesselh(k, y, prec);
     353          21 :         if (v) y = sertoser(y, t + k*v); /* lower precision */
     354          21 :         y = gdiv(a, gpowgs(y, k));
     355          21 :         T = cgetg(k+1, t_VECSMALL);
     356         168 :         for (i = 1; i <= k; i++) T[i] = 2*i+1;
     357          21 :         y = gmul(y, zv_prod_Z(T));
     358             :       }
     359          28 :       return gerepileupto(av, y);
     360             :     }
     361             :   }
     362           0 :   pari_err_TYPE("besseljh",z);
     363             :   return NULL; /* LCOV_EXCL_LINE */
     364             : }
     365             : 
     366             : static GEN
     367           0 : kbessel2(GEN nu, GEN x, long prec)
     368             : {
     369           0 :   pari_sp av = avma;
     370             :   GEN p1, x2, a;
     371             : 
     372           0 :   if (typ(x)==t_REAL) prec = realprec(x);
     373           0 :   x2 = gshift(x,1);
     374           0 :   a = gtofp(gaddgs(gshift(nu,1), 1), prec);
     375           0 :   p1 = hyperu(gshift(a,-1),a,x2,prec);
     376           0 :   p1 = gmul(gmul(p1,gpow(x2,nu,prec)), sqrtr(mppi(prec)));
     377           0 :   return gerepileupto(av, gmul(p1,gexp(gneg(x),prec)));
     378             : }
     379             : 
     380             : /* special case of hyperu */
     381             : static GEN
     382          14 : kbessel1(GEN nu, GEN gx, long prec)
     383             : {
     384             :   GEN x, y, zf, r, u, pi, nu2;
     385             :   long bit, l, k, k2, n2, n;
     386             :   pari_sp av;
     387             : 
     388          14 :   if (typ(nu)==t_COMPLEX) return kbessel2(nu,gx,prec);
     389          14 :   l = (typ(gx)==t_REAL)? realprec(gx): prec;
     390          14 :   y = cgetr(l); av = avma;
     391          14 :   x = gtofp(gx, l);
     392          14 :   nu = gtofp(nu,l); nu2 = sqrr(nu);
     393          14 :   shiftr_inplace(nu2,2); togglesign(nu2); /* nu2 = -4nu^2 */
     394          14 :   n = (long) (prec2nbits_mul(l,M_LN2) + M_PI*fabs(rtodbl(nu))) / 2;
     395          14 :   bit = prec2nbits(l) - 1;
     396          14 :   l += EXTRAPRECWORD;
     397          14 :   pi = mppi(l); n2 = n<<1; r = gmul2n(x,1);
     398          14 :   if (cmprs(x, n) < 0)
     399             :   {
     400          14 :     pari_sp av2 = avma;
     401          14 :     GEN q, v, c, s = real_1(l), t = real_0(l);
     402        1246 :     for (k = n2, k2 = 2*n2-1; k > 0; k--, k2 -= 2)
     403             :     {
     404        1232 :       GEN ak = divri(addri(nu2, sqru(k2)), mulss(n2<<2, -k));
     405        1232 :       s = addsr(1, mulrr(ak,s));
     406        1232 :       t = addsr(k2,mulrr(ak,t));
     407        1232 :       if (gc_needed(av2,3)) gerepileall(av2, 2, &s,&t);
     408             :     }
     409          14 :     shiftr_inplace(t, -1);
     410          14 :     q = utor(n2, l);
     411          14 :     zf = sqrtr(divru(pi,n2));
     412          14 :     u = gprec_wensure(mulrr(zf, s), l);
     413          14 :     v = gprec_wensure(divrs(addrr(mulrr(t,zf),mulrr(u,nu)),-n2), l);
     414             :     for(;;)
     415         301 :     {
     416         315 :       GEN p1, e, f, d = real_1(l);
     417             :       pari_sp av3;
     418         315 :       c = divur(5,q); if (expo(c) >= -1) c = real2n(-1,l);
     419         315 :       p1 = subsr(1, divrr(r,q)); if (cmprr(c,p1)>0) c = p1;
     420         315 :       togglesign(c); av3 = avma;
     421         315 :       e = u;
     422         315 :       f = v;
     423         315 :       for (k = 1;; k++)
     424       34475 :       {
     425       34790 :         GEN w = addrr(gmul2n(mulur(2*k-1,u), -1), mulrr(subrs(q,k),v));
     426       34790 :         w = addrr(w, mulrr(nu, subrr(u,gmul2n(v,1))));
     427       34790 :         u = divru(mulrr(q,v), k);
     428       34790 :         v = divru(w,k);
     429       34790 :         d = mulrr(d,c);
     430       34790 :         e = addrr(e, mulrr(d,u));
     431       34790 :         f = addrr(f, p1 = mulrr(d,v));
     432       34790 :         if (expo(p1) - expo(f) <= 1-prec2nbits(realprec(p1))) break;
     433       34475 :         if (gc_needed(av3,3)) gerepileall(av3,5,&u,&v,&d,&e,&f);
     434             :       }
     435         315 :       u = e;
     436         315 :       v = f;
     437         315 :       q = mulrr(q, addrs(c,1));
     438         315 :       if (expo(r) - expo(subrr(q,r)) >= bit) break;
     439         301 :       gerepileall(av2, 3, &u,&v,&q);
     440             :     }
     441          14 :     u = mulrr(u, gpow(divru(x,n),nu,prec));
     442             :   }
     443             :   else
     444             :   {
     445           0 :     GEN s, zz = ginv(gmul2n(r,2));
     446           0 :     pari_sp av2 = avma;
     447           0 :     s = real_1(l);
     448           0 :     for (k = n2, k2 = 2*n2-1; k > 0; k--, k2 -= 2)
     449             :     {
     450           0 :       GEN ak = divru(mulrr(addri(nu2, sqru(k2)), zz), k);
     451           0 :       s = subsr(1, mulrr(ak,s));
     452           0 :       if (gc_needed(av2,3)) s = gerepileuptoleaf(av2, s);
     453             :     }
     454           0 :     zf = sqrtr(divrr(pi,r));
     455           0 :     u = mulrr(s, zf);
     456             :   }
     457          14 :   affrr(mulrr(u, mpexp(negr(x))), y);
     458          14 :   set_avma(av); return y;
     459             : }
     460             : 
     461             : /*   sum_{k=0}^m Z^k (H(k)+H(k+n)) / (k! (k+n)!)
     462             :  * + sum_{k=0}^{n-1} (-Z)^(k-n) (n-k-1)!/k!   with Z := (-1)^flag*z^2/4.
     463             :  * Warning: contrary to _jbessel, no n! in front.
     464             :  * When flag > 1, compute exactly the H(k) and factorials (slow) */
     465             : static GEN
     466        2800 : _kbessel1(long n, GEN z, long flag, long m, long prec)
     467             : {
     468             :   GEN Z, p1, p2, s, H;
     469             :   pari_sp av;
     470             :   long k;
     471             : 
     472        2800 :   Z = gmul2n(gsqr(z),-2); if (flag & 1) Z = gneg(Z);
     473        2800 :   if (typ(z) == t_SER)
     474             :   {
     475          98 :     long v = valp(z);
     476          98 :     if (v < 0) pari_err_DOMAIN("_kbessel1","valuation", "<", gen_0, z);
     477          84 :     k = lg(Z)-2 - v;
     478          84 :     if (v == 0) pari_err_IMPL("Bessel K around a!=0");
     479          84 :     if (k <= 0) return gadd(gen_1, zeroser(varn(z), 2*v));
     480          84 :     setlg(Z, k+2);
     481             :   }
     482        2786 :   H = cgetg(m+n+2,t_VEC); gel(H,1) = gen_0;
     483        2786 :   if (flag <= 1)
     484             :   {
     485        2702 :     gel(H,2) = s = real_1(prec);
     486       98840 :     for (k=2; k<=m+n; k++) gel(H,k+1) = s = divru(addsr(1,mulur(k,s)),k);
     487             :   }
     488             :   else
     489             :   {
     490          84 :     gel(H,2) = s = gen_1;
     491         812 :     for (k=2; k<=m+n; k++) gel(H,k+1) = s = gdivgs(gaddsg(1,gmulsg(k,s)),k);
     492             :   }
     493        2786 :   s = gadd(gel(H,m+1), gel(H,m+n+1));
     494        2786 :   av = avma;
     495      101066 :   for (k=m; k>0; k--)
     496             :   {
     497       98280 :     s = gadd(gadd(gel(H,k),gel(H,k+n)),gdiv(gmul(Z,s),mulss(k,k+n)));
     498       98280 :     if (gc_needed(av,1))
     499             :     {
     500           0 :       if (DEBUGMEM>1) pari_warn(warnmem,"_kbessel1");
     501           0 :       s = gerepileupto(av, s);
     502             :     }
     503             :   }
     504        2786 :   p1 = (flag <= 1) ? mpfactr(n,prec) : mpfact(n);
     505        2786 :   s = gdiv(s,p1);
     506        2786 :   if (n)
     507             :   {
     508         364 :     Z = gneg(ginv(Z));
     509         364 :     p2 = gmulsg(n, gdiv(Z,p1));
     510         364 :     s = gadd(s,p2);
     511        1372 :     for (k=n-1; k>0; k--)
     512             :     {
     513        1008 :       p2 = gmul(p2, gmul(mulss(k,n-k),Z));
     514        1008 :       s = gadd(s,p2);
     515             :     }
     516             :   }
     517        2786 :   return s;
     518             : }
     519             : 
     520             : /* flag = 0: K / flag = 1: Y */
     521             : static GEN
     522        4347 : kbesselintern(GEN n, GEN z, long flag, long prec)
     523             : {
     524        4347 :   const char *f = flag? "besseln": "besselk";
     525        4347 :   const long flK = (flag == 0);
     526             :   long i, k, ki, lim, precnew, fl2, ex, bit;
     527        4347 :   pari_sp av = avma;
     528             :   GEN p1, p2, y, p3, pp, pm, s, c;
     529             :   double az;
     530             : 
     531        4347 :   switch(typ(z))
     532             :   {
     533           0 :     case t_QUAD: z = gtofp(z, prec); /* don't bother */
     534        3969 :     case t_INT: case t_FRAC: case t_REAL: case t_COMPLEX:
     535        3969 :       if (gequal0(z)) pari_err_DOMAIN(f, "argument", "=", gen_0, z);
     536        3969 :       i = precision(z); if (i) prec = i;
     537        3969 :       i = precision(n); if (i && prec > i) prec = i;
     538        3969 :       bit = prec2nbits(prec);
     539        3969 :       if (bessel_asymp(z, bit))
     540             :       {
     541         924 :         GEN R = flK? besskasymp(n, z, bit): bessyasymp(n, z, bit);
     542         924 :         if (typ(R) == t_COMPLEX && isexactzero(imag_i(n))
     543         217 :                                 && gsigne(real_i(z)) > 0
     544          77 :                                 && isexactzero(imag_i(z))) R = gcopy(gel(R,1));
     545         924 :         return gerepileupto(av, R);
     546             :       }
     547             :       /* heuristic threshold */
     548        3045 :       if (flK && !gequal0(n) && gexpo(z) > bit/16 + gexpo(n))
     549          14 :         return kbessel1(n,z,prec);
     550        3024 :       az = dblmodulus(z); precnew = prec;
     551        3024 :       if (az >= 1) precnew += 1 + nbits2extraprec((long)((flK?2*az:az)/M_LN2));
     552        3024 :       z = gtofp(z, precnew);
     553        3024 :       if (issmall(n,&ki))
     554             :       {
     555        2702 :         GEN z2 = gmul2n(z, -1);
     556        2702 :         double B, L = HALF_E * az;
     557        2702 :         k = labs(ki);
     558        2702 :         B = prec2nbits_mul(prec,M_LN2/2) / L;
     559        2702 :         if (flK) B += 0.367879; /* exp(-1) */
     560        2702 :         lim = bessel_get_lim(B, L);
     561        2702 :         p1 = gmul(gpowgs(z2,k), _kbessel1(k,z,flag,lim,precnew));
     562        2702 :         p2 = gadd(mpeuler(precnew), glog(z2,precnew));
     563        2702 :         p3 = jbesselintern(stoi(k),z,flag,precnew);
     564        2702 :         p2 = gsub(gmul2n(p1,-1),gmul(p2,p3));
     565        2702 :         p2 = gprec_wtrunc(p2, prec);
     566        2702 :         if (flK) {
     567        2408 :           if (k & 1) p2 = gneg(p2);
     568             :         }
     569             :         else
     570             :         {
     571         294 :           p2 = gdiv(p2, Pi2n(-1,prec));
     572         294 :           if (ki >= 0 || (k&1)==0) p2 = gneg(p2);
     573             :         }
     574        2702 :         return gerepilecopy(av, p2);
     575             :       }
     576             : 
     577         259 :       n = gtofp(n, precnew);
     578         259 :       gsincos(gmul(n,mppi(precnew)), &s,&c,precnew);
     579         259 :       ex = gexpo(s);
     580         259 :       if (ex < 0) precnew += nbits2extraprec(flK? -2*ex: -ex);
     581         259 :       if (i && i < precnew) {
     582          84 :         n = gtofp(n,precnew);
     583          84 :         z = gtofp(z,precnew);
     584          84 :         gsincos(gmul(n,mppi(precnew)), &s,&c,precnew);
     585             :       }
     586             : 
     587         259 :       pp = jbesselintern(n,      z,flag,precnew);
     588         259 :       pm = jbesselintern(gneg(n),z,flag,precnew);
     589         259 :       if (flK)
     590          70 :         p1 = gmul(gsub(pm,pp), Pi2n(-1,precnew));
     591             :       else
     592         189 :         p1 = gsub(gmul(c,pp),pm);
     593         259 :       p1 = gdiv(p1, s);
     594         259 :       return gerepilecopy(av, gprec_wtrunc(p1,prec));
     595             : 
     596         112 :     case t_VEC: case t_COL: case t_MAT:
     597         112 :       return kbesselvec(n,z,flag,prec);
     598             : 
     599          28 :     case t_POLMOD:
     600          28 :       y = kbesselvec(n,polmod_to_embed(z,prec),flag,prec);
     601          28 :       return gerepileupto(av, y);
     602             : 
     603          14 :     case t_PADIC: pari_err_IMPL(stack_strcat("p-adic ",f));
     604         224 :     default:
     605         224 :       if (!(y = toser_i(z))) break;
     606         210 :       if (issmall(n,&ki))
     607             :       {
     608          98 :         k = labs(ki);
     609          98 :         return gerepilecopy(av, _kbessel1(k,y,flag+2,lg(y)-2,prec));
     610             :       }
     611          98 :       if (!issmall(gmul2n(n,1),&ki))
     612          70 :         pari_err_DOMAIN(f, "2n mod Z", "!=", gen_0,n);
     613          28 :       k = labs(ki); n = gmul2n(stoi(k),-1);
     614          28 :       fl2 = (k&3)==1;
     615          28 :       pm = jbesselintern(gneg(n),y,flag,prec);
     616          28 :       if (flK)
     617             :       {
     618           7 :         pp = jbesselintern(n,y,flag,prec);
     619           7 :         p2 = gpowgs(y,-k); if (fl2 == 0) p2 = gneg(p2);
     620           7 :         p3 = gmul2n(diviiexact(mpfact(k + 1),mpfact((k + 1) >> 1)),-(k + 1));
     621           7 :         p3 = gdivgs(gmul2n(gsqr(p3),1),k);
     622           7 :         p2 = gmul(p2,p3);
     623           7 :         p1 = gsub(pp,gmul(p2,pm));
     624             :       }
     625          21 :       else p1 = pm;
     626          28 :       return gerepileupto(av, fl2? gneg(p1): gcopy(p1));
     627             :   }
     628          14 :   pari_err_TYPE(f,z);
     629             :   return NULL; /* LCOV_EXCL_LINE */
     630             : }
     631             : 
     632             : GEN
     633        3059 : kbessel(GEN n, GEN z, long prec) { return kbesselintern(n,z,0,prec); }
     634             : GEN
     635        1120 : nbessel(GEN n, GEN z, long prec) { return kbesselintern(n,z,1,prec); }
     636             : /* J + iN */
     637             : GEN
     638         224 : hbessel1(GEN n, GEN z, long prec)
     639             : {
     640         224 :   pari_sp av = avma;
     641         224 :   GEN J = jbessel(n,z,prec);
     642         196 :   GEN N = nbessel(n,z,prec);
     643         182 :   return gerepileupto(av, gadd(J, mulcxI(N)));
     644             : }
     645             : /* J - iN */
     646             : GEN
     647         224 : hbessel2(GEN n, GEN z, long prec)
     648             : {
     649         224 :   pari_sp av = avma;
     650         224 :   GEN J = jbessel(n,z,prec);
     651         196 :   GEN N = nbessel(n,z,prec);
     652         182 :   return gerepileupto(av, gadd(J, mulcxmI(N)));
     653             : }
     654             : 
     655             : /***********************************************************************/
     656             : /**                    INCOMPLETE GAMMA FUNCTION                      **/
     657             : /***********************************************************************/
     658             : /* mx ~ |x|, b = bit accuracy */
     659             : static int
     660       16366 : gamma_use_asymp(GEN x, long b)
     661             : {
     662             :   long e;
     663       16366 :   if (is_real_t(typ(x)))
     664             :   {
     665       12348 :     pari_sp av = avma;
     666       12348 :     return gc_int(av, gcmpgs(R_abs_shallow(x), 3*b / 4) >= 0);
     667             :   }
     668        4018 :   e = gexpo(x); return e >= b || dblmodulus(x) >= 3*b / 4;
     669             : }
     670             : /* x a t_REAL */
     671             : static GEN
     672          28 : eint1r_asymp(GEN x, GEN expx, long prec)
     673             : {
     674          28 :   pari_sp av = avma, av2;
     675             :   GEN S, q, z, ix;
     676          28 :   long oldeq = LONG_MAX, esx = -prec2nbits(prec), j;
     677             : 
     678          28 :   if (realprec(x) < prec + EXTRAPREC) x = rtor(x, prec+EXTRAPREC);
     679          28 :   ix = invr(x); q = z = negr(ix);
     680          28 :   av2 = avma; S = addrs(q, 1);
     681          28 :   for (j = 2;; j++)
     682        1211 :   {
     683        1239 :     long eq = expo(q); if (eq < esx) break;
     684        1211 :     if ((j & 3) == 0)
     685             :     { /* guard against divergence */
     686         294 :       if (eq > oldeq) return gc_NULL(av); /* regressing, abort */
     687         294 :       oldeq = eq;
     688             :     }
     689        1211 :     q = mulrr(q, mulru(z, j)); S = addrr(S, q);
     690        1211 :     if (gc_needed(av2, 1)) gerepileall(av2, 2, &S, &q);
     691             :   }
     692          28 :   if (DEBUGLEVEL > 2) err_printf("eint1: using asymp\n");
     693          28 :   S = expx? divrr(S, expx): mulrr(S, mpexp(negr(x)));
     694          28 :   return gerepileuptoleaf(av, mulrr(S, ix));
     695             : }
     696             : /* cf incgam_asymp(0, x); z = -1/x
     697             :  *   exp(-x)/x * (1 + z + 2! z^2 + ...) */
     698             : static GEN
     699         105 : eint1_asymp(GEN x, GEN expx, long prec)
     700             : {
     701         105 :   pari_sp av = avma, av2;
     702             :   GEN S, q, z, ix;
     703         105 :   long oldeq = LONG_MAX, esx = -prec2nbits(prec), j;
     704             : 
     705         105 :   if (typ(x) != t_REAL) x = gtofp(x, prec+EXTRAPREC);
     706         105 :   if (typ(x) == t_REAL) return eint1r_asymp(x, expx, prec);
     707         105 :   ix = ginv(x); q = z = gneg_i(ix);
     708         105 :   av2 = avma; S = gaddgs(q, 1);
     709         105 :   for (j = 2;; j++)
     710        5824 :   {
     711        5929 :     long eq = gexpo(q); if (eq < esx) break;
     712        5824 :     if ((j & 3) == 0)
     713             :     { /* guard against divergence */
     714        1442 :       if (eq > oldeq) return gc_NULL(av); /* regressing, abort */
     715        1442 :       oldeq = eq;
     716             :     }
     717        5824 :     q = gmul(q, gmulgs(z, j)); S = gadd(S, q);
     718        5824 :     if (gc_needed(av2, 1)) gerepileall(av2, 2, &S, &q);
     719             :   }
     720         105 :   if (DEBUGLEVEL > 2) err_printf("eint1: using asymp\n");
     721         105 :   S = expx? gdiv(S, expx): gmul(S, gexp(gneg_i(x), prec));
     722         105 :   return gerepileupto(av, gmul(S, ix));
     723             : }
     724             : 
     725             : /* eint1(x) = incgam(0, x); typ(x) = t_REAL, x > 0 */
     726             : static GEN
     727        6146 : eint1p(GEN x, GEN expx)
     728             : {
     729             :   pari_sp av;
     730        6146 :   long l = realprec(x), bit = prec2nbits(l), prec, i;
     731             :   double mx;
     732             :   GEN z, S, t, H, run;
     733             : 
     734        6146 :   if (gamma_use_asymp(x, bit)
     735          28 :       && (z = eint1r_asymp(x, expx, l))) return z;
     736        6118 :   mx = rtodbl(x);
     737        6118 :   prec = l + nbits2extraprec((mx+log(mx))/M_LN2 + 10);
     738        6118 :   bit = prec2nbits(prec);
     739        6118 :   run = real_1(prec); x = rtor(x, prec);
     740        6118 :   av = avma; S = z = t = H = run;
     741      622049 :   for (i = 2; expo(S) - expo(t) <= bit; i++)
     742             :   {
     743      615931 :     H = addrr(H, divru(run,i)); /* H = sum_{k<=i} 1/k */
     744      615931 :     z = divru(mulrr(x,z), i);   /* z = x^(i-1)/i! */
     745      615931 :     t = mulrr(z, H); S = addrr(S, t);
     746      615931 :     if ((i & 0x1ff) == 0) gerepileall(av, 4, &z,&t,&S,&H);
     747             :   }
     748        6118 :   return subrr(mulrr(x, divrr(S,expx? expx: mpexp(x))),
     749             :                addrr(mplog(x), mpeuler(prec)));
     750             : }
     751             : /* eint1(x) = incgam(0, x); typ(x) = t_REAL, x < 0
     752             :  * rewritten from code contributed by Manfred Radimersky */
     753             : static GEN
     754         140 : eint1m(GEN x, GEN expx)
     755             : {
     756         140 :   GEN p1, q, S, y, z = cgetg(3, t_COMPLEX);
     757         140 :   long l  = realprec(x), n  = prec2nbits(l), j;
     758         140 :   pari_sp av = avma;
     759             : 
     760         140 :   y  = rtor(x, l + EXTRAPREC); setsigne(y,1); /* |x| */
     761         140 :   if (gamma_use_asymp(y, n))
     762             :   { /* ~eint1_asymp: asymptotic expansion */
     763          14 :     p1 = q = invr(y); S = addrs(q, 1);
     764         560 :     for (j = 2; expo(q) >= -n; j++) {
     765         546 :       q = mulrr(q, mulru(p1, j));
     766         546 :       S = addrr(S, q);
     767             :     }
     768          14 :     y  = mulrr(p1, expx? divrr(S, expx): mulrr(S, mpexp(y)));
     769             :   }
     770             :   else
     771             :   {
     772         126 :     p1 = q = S = y;
     773       24248 :     for (j = 2; expo(q) - expo(S) >= -n; j++) {
     774       24122 :       p1 = mulrr(y, divru(p1, j)); /* (-x)^j/j! */
     775       24122 :       q = divru(p1, j);
     776       24122 :       S = addrr(S, q);
     777             :     }
     778         126 :     y  = addrr(S, addrr(logr_abs(x), mpeuler(l)));
     779             :   }
     780         140 :   y = gerepileuptoleaf(av, y); togglesign(y);
     781         140 :   gel(z, 1) = y;
     782         140 :   y = mppi(l); setsigne(y, -1);
     783         140 :   gel(z, 2) = y; return z;
     784             : }
     785             : 
     786             : /* real(z*log(z)-z), z = x+iy */
     787             : static double
     788        7672 : mygamma(double x, double y)
     789             : {
     790        7672 :   if (x == 0.) return -(M_PI/2)*fabs(y);
     791        7672 :   return (x/2)*log(x*x+y*y)-x-y*atan(y/x);
     792             : }
     793             : 
     794             : /* x^s exp(-x) */
     795             : static GEN
     796       10010 : expmx_xs(GEN s, GEN x, GEN logx, long prec)
     797             : {
     798             :   GEN z;
     799       10010 :   long ts = typ(s);
     800       10010 :   if (ts == t_INT || (ts == t_FRAC && absequaliu(gel(s,2), 2)))
     801        5243 :     z = gmul(gexp(gneg(x), prec), gpow(x, s, prec));
     802             :   else
     803        4767 :     z = gexp(gsub(gmul(s, logx? logx: glog(x,prec+EXTRAPREC)), x), prec);
     804       10010 :   return z;
     805             : }
     806             : 
     807             : /* Not yet: doesn't work at low accuracy
     808             : #define INCGAM_CF
     809             : */
     810             : 
     811             : #ifdef INCGAM_CF
     812             : /* Is s very close to a non-positive integer ? */
     813             : static int
     814             : isgammapole(GEN s, long bitprec)
     815             : {
     816             :   pari_sp av = avma;
     817             :   GEN t = imag_i(s);
     818             :   long e, b = bitprec - 10;
     819             : 
     820             :   if (gexpo(t) > - b) return 0;
     821             :   s = real_i(s);
     822             :   if (gsigne(s) > 0 && gexpo(s) > -b) return 0;
     823             :   (void)grndtoi(s, &e); return gc_bool(av, e < -b);
     824             : }
     825             : 
     826             : /* incgam using the continued fraction. x a t_REAL or t_COMPLEX, mx ~ |x|.
     827             :  * Assume precision(s), precision(x) >= prec */
     828             : static GEN
     829             : incgam_cf(GEN s, GEN x, double mx, long prec)
     830             : {
     831             :   GEN ms, y, S;
     832             :   long n, i, j, LS, bitprec = prec2nbits(prec);
     833             :   double rs, is, m;
     834             : 
     835             :   if (typ(s) == t_COMPLEX)
     836             :   {
     837             :     rs = gtodouble(gel(s,1));
     838             :     is = gtodouble(gel(s,2));
     839             :   }
     840             :   else
     841             :   {
     842             :     rs = gtodouble(s);
     843             :     is = 0.;
     844             :   }
     845             :   if (isgammapole(s, bitprec)) LS = 0;
     846             :   else
     847             :   {
     848             :     double bit,  LGS = mygamma(rs,is);
     849             :     LS = LGS <= 0 ? 0: ceil(LGS);
     850             :     bit = (LGS - (rs-1)*log(mx) + mx)/M_LN2;
     851             :     if (bit > 0)
     852             :     {
     853             :       prec += nbits2extraprec((long)bit);
     854             :       x = gtofp(x, prec);
     855             :       if (isinexactreal(s)) s = gtofp(s, prec);
     856             :     }
     857             :   }
     858             :   /* |ln(2*gamma(s)*sin(s*Pi))| <= ln(2) + |lngamma(s)| + |Im(s)*Pi|*/
     859             :   m = bitprec*M_LN2 + LS + M_LN2 + fabs(is)*M_PI + mx;
     860             :   if (rs < 1) m += (1 - rs)*log(mx);
     861             :   m /= 4;
     862             :   n = (long)(1 + m*m/mx);
     863             :   y = expmx_xs(gsubgs(s,1), x, NULL, prec);
     864             :   if (rs >= 0 && bitprec >= 512)
     865             :   {
     866             :     GEN A = cgetg(n+1, t_VEC), B = cgetg(n+1, t_VEC);
     867             :     ms = gsubsg(1, s);
     868             :     for (j = 1; j <= n; ++j)
     869             :     {
     870             :       gel(A,j) = ms;
     871             :       gel(B,j) = gmulsg(j, gsubgs(s,j));
     872             :       ms = gaddgs(ms, 2);
     873             :     }
     874             :     S = contfraceval_inv(mkvec2(A,B), x, -1);
     875             :   }
     876             :   else
     877             :   {
     878             :     GEN x_s = gsub(x, s);
     879             :     pari_sp av2 = avma;
     880             :     S = gdiv(gsubgs(s,n), gaddgs(x_s,n<<1));
     881             :     for (i=n-1; i >= 1; i--)
     882             :     {
     883             :       S = gdiv(gsubgs(s,i), gadd(gaddgs(x_s,i<<1),gmulsg(i,S)));
     884             :       if (gc_needed(av2,3))
     885             :       {
     886             :         if(DEBUGMEM>1) pari_warn(warnmem,"incgam_cf");
     887             :         S = gerepileupto(av2, S);
     888             :       }
     889             :     }
     890             :     S = gaddgs(S,1);
     891             :   }
     892             :   return gmul(y, S);
     893             : }
     894             : #endif
     895             : 
     896             : static double
     897        5586 : findextraincgam(GEN s, GEN x)
     898             : {
     899        5586 :   double sig = gtodouble(real_i(s)), t = gtodouble(imag_i(s));
     900        5586 :   double xr = gtodouble(real_i(x)), xi = gtodouble(imag_i(x));
     901        5586 :   double exd = 0., Nx = xr*xr + xi*xi, D = Nx - t*t;
     902             :   long n;
     903             : 
     904        5586 :   if (xr < 0)
     905             :   {
     906         819 :     long ex = gexpo(x);
     907         819 :     if (ex > 0 && ex > gexpo(s)) exd = sqrt(Nx)*log(Nx)/2; /* |x| log |x| */
     908             :   }
     909        5586 :   if (D <= 0.) return exd;
     910        4144 :   n = (long)(sqrt(D)-sig);
     911        4144 :   if (n <= 0) return exd;
     912        1491 :   return maxdd(exd, (n*log(Nx)/2 - mygamma(sig+n, t) + mygamma(sig, t)) / M_LN2);
     913             : }
     914             : 
     915             : /* use exp(-x) * (x^s/s) * sum_{k >= 0} x^k / prod(i=1, k, s+i) */
     916             : static GEN
     917        5593 : incgamc_i(GEN s, GEN x, long *ptexd, long prec)
     918             : {
     919             :   GEN S, t, y;
     920             :   long l, n, i, exd;
     921        5593 :   pari_sp av = avma, av2;
     922             : 
     923        5593 :   if (gequal0(x))
     924             :   {
     925           7 :     if (ptexd) *ptexd = 0.;
     926           7 :     return gtofp(x, prec);
     927             :   }
     928        5586 :   l = precision(x);
     929        5586 :   if (!l) l = prec;
     930        5586 :   n = -prec2nbits(l)-1;
     931        5586 :   exd = (long)findextraincgam(s, x);
     932        5586 :   if (ptexd) *ptexd = exd;
     933        5586 :   if (exd > 0)
     934             :   {
     935        1372 :     long p = l + nbits2extraprec(exd);
     936        1372 :     x = gtofp(x, p);
     937        1372 :     if (isinexactreal(s)) s = gtofp(s, p);
     938             :   }
     939        4214 :   else x = gtofp(x, l+EXTRAPREC);
     940        5586 :   av2 = avma;
     941        5586 :   S = gdiv(x, gaddsg(1,s));
     942        5586 :   t = gaddsg(1, S);
     943      751764 :   for (i=2; gexpo(S) >= n; i++)
     944             :   {
     945      746178 :     S = gdiv(gmul(x,S), gaddsg(i,s)); /* x^i / ((s+1)...(s+i)) */
     946      746178 :     t = gadd(S,t);
     947      746178 :     if (gc_needed(av2,3))
     948             :     {
     949           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"incgamc");
     950           0 :       gerepileall(av2, 2, &S, &t);
     951             :     }
     952             :   }
     953        5586 :   y = expmx_xs(s, x, NULL, prec);
     954        5586 :   return gerepileupto(av, gmul(gdiv(y,s), t));
     955             : }
     956             : 
     957             : GEN
     958        1414 : incgamc(GEN s, GEN x, long prec)
     959        1414 : { return incgamc_i(s, x, NULL, prec); }
     960             : 
     961             : /* incgamma using asymptotic expansion:
     962             :  *   exp(-x)x^(s-1)(1 + (s-1)/x + (s-1)(s-2)/x^2 + ...) */
     963             : static GEN
     964        2716 : incgam_asymp(GEN s, GEN x, long prec)
     965             : {
     966        2716 :   pari_sp av = avma, av2;
     967             :   GEN S, q, cox, invx;
     968        2716 :   long oldeq = LONG_MAX, eq, esx, j;
     969        2716 :   int flint = (typ(s) == t_INT && signe(s) > 0);
     970             : 
     971        2716 :   x = gtofp(x,prec+EXTRAPREC);
     972        2716 :   invx = ginv(x);
     973        2716 :   esx = -prec2nbits(prec);
     974        2716 :   av2 = avma;
     975        2716 :   q = gmul(gsubgs(s,1), invx);
     976        2716 :   S = gaddgs(q, 1);
     977        2716 :   for (j = 2;; j++)
     978             :   {
     979      123879 :     eq = gexpo(q); if (eq < esx) break;
     980      121282 :     if (!flint && (j & 3) == 0)
     981             :     { /* guard against divergence */
     982       15778 :       if (eq > oldeq) return gc_NULL(av); /* regressing, abort */
     983       15659 :       oldeq = eq;
     984             :     }
     985      121163 :     q = gmul(q, gmul(gsubgs(s,j), invx));
     986      121163 :     S = gadd(S, q);
     987      121163 :     if (gc_needed(av2, 1)) gerepileall(av2, 2, &S, &q);
     988             :   }
     989        2597 :   if (DEBUGLEVEL > 2) err_printf("incgam: using asymp\n");
     990        2597 :   cox = expmx_xs(gsubgs(s,1), x, NULL, prec);
     991        2597 :   return gerepileupto(av, gmul(cox, S));
     992             : }
     993             : 
     994             : /* gasx = incgam(s-n,x). Compute incgam(s,x)
     995             :  * = (s-1)(s-2)...(s-n)gasx + exp(-x)x^(s-1) *
     996             :  *   (1 + (s-1)/x + ... + (s-1)(s-2)...(s-n+1)/x^(n-1)) */
     997             : static GEN
     998         546 : incgam_asymp_partial(GEN s, GEN x, GEN gasx, long n, long prec)
     999             : {
    1000             :   pari_sp av;
    1001         546 :   GEN S, q, cox, invx, s1 = gsubgs(s, 1), sprod;
    1002             :   long j;
    1003         546 :   cox = expmx_xs(s1, x, NULL, prec);
    1004         546 :   if (n == 1) return gadd(cox, gmul(s1, gasx));
    1005         546 :   invx = ginv(x);
    1006         546 :   av = avma;
    1007         546 :   q = gmul(s1, invx);
    1008         546 :   S = gaddgs(q, 1);
    1009       52164 :   for (j = 2; j < n; j++)
    1010             :   {
    1011       51618 :     q = gmul(q, gmul(gsubgs(s, j), invx));
    1012       51618 :     S = gadd(S, q);
    1013       51618 :     if (gc_needed(av, 2)) gerepileall(av, 2, &S, &q);
    1014             :   }
    1015         546 :   sprod = gmul(gmul(q, gpowgs(x, n-1)), gsubgs(s, n));
    1016         546 :   return gadd(gmul(cox, S), gmul(sprod, gasx));
    1017             : }
    1018             : 
    1019             : /* Assume s != 0; called when Re(s) <= 1/2 */
    1020             : static GEN
    1021        2401 : incgamspec(GEN s, GEN x, GEN g, long prec)
    1022             : {
    1023        2401 :   GEN q, S, cox = gen_0, P, sk, S1, S2, S3, F3, logx, mx;
    1024        2401 :   long n, esk, E, k = itos(ground(gneg(real_i(s)))); /* >= 0 */
    1025             : 
    1026        2401 :   if (k && gexpo(x) > 0)
    1027             :   {
    1028         245 :     GEN xk = gdivgs(x, k);
    1029         245 :     long bitprec = prec2nbits(prec);
    1030         245 :     double d = (gexpo(xk) > bitprec)? bitprec*M_LN2: log(dblmodulus(xk));
    1031         245 :     d = k * (d + 1) / M_LN2;
    1032         245 :     if (d > 0) prec += nbits2extraprec((long)d);
    1033         245 :     if (isinexactreal(s)) s = gtofp(s, prec);
    1034             :   }
    1035        2401 :   x = gtofp(x, maxss(precision(x), prec) + EXTRAPREC);
    1036        2401 :   sk = gaddgs(s, k); /* |Re(sk)| <= 1/2 */
    1037        2401 :   logx = glog(x, prec);
    1038        2401 :   mx = gneg(x);
    1039        2401 :   if (k == 0) { S = gen_0; P = gen_1; }
    1040             :   else
    1041             :   {
    1042             :     long j;
    1043         854 :     q = ginv(s); S = q; P = s;
    1044       16926 :     for (j = 1; j < k; j++)
    1045             :     {
    1046       16072 :       GEN sj = gaddgs(s, j);
    1047       16072 :       q = gmul(q, gdiv(x, sj));
    1048       16072 :       S = gadd(S, q);
    1049       16072 :       P = gmul(P, sj);
    1050             :     }
    1051         854 :     cox = expmx_xs(s, x, logx, prec); /* x^s exp(-x) */
    1052         854 :     S = gmul(S, gneg(cox));
    1053             :   }
    1054        2401 :   if (k && gequal0(sk))
    1055         175 :     return gadd(S, gdiv(eint1(x, prec), P));
    1056        2226 :   esk = gexpo(sk);
    1057        2226 :   if (esk > -7)
    1058             :   {
    1059        1015 :     GEN a, b, PG = gmul(sk, P);
    1060        1015 :     if (g) g = gmul(g, PG);
    1061        1015 :     a = incgam0(gaddgs(sk,1), x, g, prec);
    1062        1015 :     if (k == 0) cox = expmx_xs(s, x, logx, prec);
    1063        1015 :     b = gmul(gpowgs(x, k), cox);
    1064        1015 :     return gadd(S, gdiv(gsub(a, b), PG));
    1065             :   }
    1066        1211 :   E = prec2nbits(prec) + 1;
    1067        1211 :   if (gexpo(x) > 0)
    1068             :   {
    1069         420 :     long X = (long)(dblmodulus(x)/M_LN2);
    1070         420 :     prec += 2*nbits2extraprec(X);
    1071         420 :     x = gtofp(x, prec); mx = gneg(x);
    1072         420 :     logx = glog(x, prec); sk = gtofp(sk, prec);
    1073         420 :     E += X;
    1074             :   }
    1075        1211 :   if (isinexactreal(sk)) sk = gtofp(sk, prec+EXTRAPREC);
    1076             :   /* |sk| < 2^-7 is small, guard against cancellation */
    1077        1211 :   F3 = gexpm1(gmul(sk, logx), prec);
    1078             :   /* ( gamma(1+sk) - exp(sk log(x))) ) / sk */
    1079        1211 :   S1 = gdiv(gsub(ggamma1m1(sk, prec+EXTRAPREC), F3), sk);
    1080        1211 :   q = x; S3 = gdiv(x, gaddsg(1,sk));
    1081      255523 :   for (n = 2; gexpo(q) - gexpo(S3) > -E; ++n)
    1082             :   {
    1083      254312 :     q = gmul(q, gdivgs(mx, n));
    1084      254312 :     S3 = gadd(S3, gdiv(q, gaddsg(n, sk)));
    1085             :   }
    1086        1211 :   S2 = gadd(gadd(S1, S3), gmul(F3, S3));
    1087        1211 :   return gadd(S, gdiv(S2, P));
    1088             : }
    1089             : 
    1090             : /* return |x| */
    1091             : double
    1092     6702523 : dblmodulus(GEN x)
    1093             : {
    1094     6702523 :   if (typ(x) == t_COMPLEX)
    1095             :   {
    1096      715897 :     double a = gtodouble(gel(x,1));
    1097      715897 :     double b = gtodouble(gel(x,2));
    1098      715897 :     return sqrt(a*a + b*b);
    1099             :   }
    1100             :   else
    1101     5986626 :     return fabs(gtodouble(x));
    1102             : }
    1103             : 
    1104             : /* Driver routine. If g != NULL, assume that g=gamma(s,prec). */
    1105             : GEN
    1106       11529 : incgam0(GEN s, GEN x, GEN g, long prec)
    1107             : {
    1108             :   pari_sp av;
    1109             :   long E, l;
    1110             :   GEN z, rs, is;
    1111             : 
    1112       11529 :   if (gequal0(x)) return g? gcopy(g): ggamma(s,prec);
    1113       11529 :   if (gequal0(s)) return eint1(x, prec);
    1114        9723 :   l = precision(s); if (!l) l = prec;
    1115        9723 :   E = prec2nbits(l);
    1116        9723 :   if (gamma_use_asymp(x, E) ||
    1117        8302 :       (typ(s) == t_INT && signe(s) > 0 && gexpo(x) >= expi(s)))
    1118        2716 :     if ((z = incgam_asymp(s, x, l))) return z;
    1119        7126 :   av = avma; E++;
    1120        7126 :   rs = real_i(s);
    1121        7126 :   is = imag_i(s);
    1122             : #ifdef INCGAM_CF
    1123             :   /* Can one use continued fraction ? */
    1124             :   if (gequal0(is) && gequal0(imag_i(x)) && gsigne(x) > 0)
    1125             :   {
    1126             :     double sd = gtodouble(rs), LB, UB;
    1127             :     double xd = gtodouble(real_i(x));
    1128             :     if (sd > 0) {
    1129             :       LB = 15 + 0.1205*E;
    1130             :       UB = 5 + 0.752*E;
    1131             :     } else {
    1132             :       LB = -6 + 0.1205*E;
    1133             :       UB = 5 + 0.752*E + fabs(sd)/54.;
    1134             :     }
    1135             :     if (xd >= LB && xd <= UB)
    1136             :     {
    1137             :       if (DEBUGLEVEL > 2) err_printf("incgam: using continued fraction\n");
    1138             :       return gerepileupto(av, incgam_cf(s, x, xd, prec));
    1139             :     }
    1140             :   }
    1141             : #endif
    1142        7126 :   if (gsigne(rs) > 0 && gexpo(rs) >= -1)
    1143             :   { /* use complementary incomplete gamma */
    1144        4725 :     long n, egs, exd, precg, es = gexpo(s);
    1145        4725 :     if (es < 0) {
    1146         581 :       l += nbits2extraprec(-es) + 1;
    1147         581 :       x = gtofp(x, l);
    1148         581 :       if (isinexactreal(s)) s = gtofp(s, l);
    1149             :     }
    1150        4725 :     n = itos(gceil(rs));
    1151        4725 :     if (n > 100)
    1152             :     {
    1153             :       GEN gasx;
    1154         546 :       n -= 100;
    1155         546 :       if (es > 0)
    1156             :       {
    1157         546 :         es = mygamma(gtodouble(rs) - n, gtodouble(is)) / M_LN2;
    1158         546 :         if (es > 0)
    1159             :         {
    1160         546 :           l += nbits2extraprec(es);
    1161         546 :           x = gtofp(x, l);
    1162         546 :           if (isinexactreal(s)) s = gtofp(s, l);
    1163             :         }
    1164             :       }
    1165         546 :       gasx = incgam0(gsubgs(s, n), x, NULL, prec);
    1166         546 :       return gerepileupto(av, incgam_asymp_partial(s, x, gasx, n, prec));
    1167             :     }
    1168        4179 :     if (DEBUGLEVEL > 2) err_printf("incgam: using power series 1\n");
    1169             :     /* egs ~ expo(gamma(s)) */
    1170        4179 :     precg = g? precision(g): 0;
    1171        4179 :     egs = g? gexpo(g): (long)(mygamma(gtodouble(rs), gtodouble(is)) / M_LN2);
    1172        4179 :     if (egs > 0) {
    1173        1946 :       l += nbits2extraprec(egs) + 1;
    1174        1946 :       x = gtofp(x, l);
    1175        1946 :       if (isinexactreal(s)) s = gtofp(s, l);
    1176        1946 :       if (precg < l) g = NULL;
    1177             :     }
    1178        4179 :     z = incgamc_i(s, x, &exd, l);
    1179        4179 :     if (exd > 0)
    1180             :     {
    1181         896 :       l += nbits2extraprec(exd);
    1182         896 :       if (isinexactreal(s)) s = gtofp(s, l);
    1183         896 :       if (precg < l) g = NULL;
    1184             :     }
    1185             :     else
    1186             :     { /* gamma(s) negligible ? Compute to lower accuracy */
    1187        3283 :       long e = gexpo(z) - egs;
    1188        3283 :       if (e > 3)
    1189             :       {
    1190         420 :         E -= e;
    1191         420 :         if (E <= 0) g = gen_0; else if (!g) g = ggamma(s, nbits2prec(E));
    1192             :       }
    1193             :     }
    1194             :     /* worry about possible cancellation */
    1195        4179 :     if (!g) g = ggamma(s, maxss(l,precision(z)));
    1196        4179 :     return gerepileupto(av, gsub(g,z));
    1197             :   }
    1198        2401 :   if (DEBUGLEVEL > 2) err_printf("incgam: using power series 2\n");
    1199        2401 :   return gerepileupto(av, incgamspec(s, x, g, l));
    1200             : }
    1201             : 
    1202             : GEN
    1203        1106 : incgam(GEN s, GEN x, long prec) { return incgam0(s, x, NULL, prec); }
    1204             : 
    1205             : /* x a t_REAL */
    1206             : GEN
    1207        2114 : mpeint1(GEN x, GEN expx)
    1208             : {
    1209        2114 :   long s = signe(x);
    1210             :   pari_sp av;
    1211             :   GEN z;
    1212        2114 :   if (!s) pari_err_DOMAIN("eint1", "x","=",gen_0, x);
    1213        2107 :   if (s < 0) return eint1m(x, expx);
    1214        1967 :   z = cgetr(realprec(x));
    1215        1967 :   av = avma; affrr(eint1p(x, expx), z);
    1216        1967 :   set_avma(av); return z;
    1217             : }
    1218             : 
    1219             : static GEN
    1220         357 : cxeint1(GEN x, long prec)
    1221             : {
    1222         357 :   pari_sp av = avma, av2;
    1223             :   GEN q, S, run, z, H;
    1224         357 :   long n, E = prec2nbits(prec);
    1225             : 
    1226         357 :   if (gamma_use_asymp(x, E) && (z = eint1_asymp(x, NULL, prec))) return z;
    1227         252 :   E++;
    1228         252 :   if (gexpo(x) > 0)
    1229             :   { /* take cancellation into account, log2(\sum |x|^n / n!) = |x| / log(2) */
    1230          42 :     double dbx = dblmodulus(x);
    1231          42 :     long X = (long)((dbx + log(dbx))/M_LN2 + 10);
    1232          42 :     prec += nbits2extraprec(X);
    1233          42 :     x = gtofp(x, prec); E += X;
    1234             :   }
    1235         252 :   if (DEBUGLEVEL > 2) err_printf("eint1: using power series\n");
    1236         252 :   run = real_1(prec);
    1237         252 :   av2 = avma;
    1238         252 :   S = z = q = H = run;
    1239       48384 :   for (n = 2; gexpo(q) - gexpo(S) >= -E; n++)
    1240             :   {
    1241       48132 :     H = addrr(H, divru(run, n)); /* H = sum_{k<=n} 1/k */
    1242       48132 :     z = gdivgs(gmul(x,z), n);   /* z = x^(n-1)/n! */
    1243       48132 :     q = gmul(z, H); S = gadd(S, q);
    1244       48132 :     if ((n & 0x1ff) == 0) gerepileall(av2, 4, &z, &q, &S, &H);
    1245             :   }
    1246         252 :   S = gmul(gmul(x, S), gexp(gneg_i(x), prec));
    1247         252 :   return gerepileupto(av, gsub(S, gadd(glog(x, prec), mpeuler(prec))));
    1248             : }
    1249             : 
    1250             : GEN
    1251        2471 : eint1(GEN x, long prec)
    1252             : {
    1253        2471 :   switch(typ(x))
    1254             :   {
    1255         357 :     case t_COMPLEX: return cxeint1(x, prec);
    1256        1729 :     case t_REAL: break;
    1257         385 :     default: x = gtofp(x, prec);
    1258             :   }
    1259        2114 :   return mpeint1(x,NULL);
    1260             : }
    1261             : 
    1262             : GEN
    1263          49 : veceint1(GEN C, GEN nmax, long prec)
    1264             : {
    1265          49 :   if (!nmax) return eint1(C,prec);
    1266           7 :   if (typ(nmax) != t_INT) pari_err_TYPE("veceint1",nmax);
    1267           7 :   if (typ(C) != t_REAL) {
    1268           7 :     C = gtofp(C, prec);
    1269           7 :     if (typ(C) != t_REAL) pari_err_TYPE("veceint1",C);
    1270             :   }
    1271           7 :   if (signe(C) <= 0) pari_err_DOMAIN("veceint1", "argument", "<=", gen_0,C);
    1272           7 :   return mpveceint1(C, NULL, itos(nmax));
    1273             : }
    1274             : 
    1275             : /* j > 0, a t_REAL. Return sum_{m >= 0} a^m / j(j+1)...(j+m)).
    1276             :  * Stop when expo(summand) < E; note that s(j-1) = (a s(j) + 1) / (j-1). */
    1277             : static GEN
    1278         259 : mp_sum_j(GEN a, long j, long E, long prec)
    1279             : {
    1280         259 :   pari_sp av = avma;
    1281         259 :   GEN q = divru(real_1(prec), j), s = q;
    1282             :   long m;
    1283         259 :   for (m = 0;; m++)
    1284             :   {
    1285        4871 :     if (expo(q) < E) break;
    1286        4612 :     q = mulrr(q, divru(a, m+j));
    1287        4612 :     s = addrr(s, q);
    1288             :   }
    1289         259 :   return gerepileuptoleaf(av, s);
    1290             : }
    1291             : /* Return the s_a(j), j <= J */
    1292             : static GEN
    1293         259 : sum_jall(GEN a, long J, long prec)
    1294             : {
    1295         259 :   GEN s = cgetg(J+1, t_VEC);
    1296         259 :   long j, E = -prec2nbits(prec) - 5;
    1297         259 :   gel(s, J) = mp_sum_j(a, J, E, prec);
    1298       11103 :   for (j = J-1; j; j--)
    1299       10844 :     gel(s,j) = divru(addrs(mulrr(a, gel(s,j+1)), 1), j);
    1300         259 :   return s;
    1301             : }
    1302             : 
    1303             : /* T a dense t_POL with t_REAL coeffs. Return T(n) [faster than poleval] */
    1304             : static GEN
    1305      481831 : rX_s_eval(GEN T, long n)
    1306             : {
    1307      481831 :   long i = lg(T)-1;
    1308      481831 :   GEN c = gel(T,i);
    1309    10743228 :   for (i--; i>=2; i--) c = gadd(mulrs(c,n),gel(T,i));
    1310      481831 :   return c;
    1311             : }
    1312             : 
    1313             : /* C>0 t_REAL, eC = exp(C). Return eint1(n*C) for 1<=n<=N. Absolute accuracy */
    1314             : GEN
    1315         266 : mpveceint1(GEN C, GEN eC, long N)
    1316             : {
    1317         266 :   const long prec = realprec(C);
    1318         266 :   long Nmin = 15; /* >= 1. E.g. between 10 and 30, but little effect */
    1319         266 :   GEN en, v, w = cgetg(N+1, t_VEC);
    1320             :   pari_sp av0;
    1321             :   double DL;
    1322             :   long n, j, jmax, jmin;
    1323         266 :   if (!N) return w;
    1324      486017 :   for (n = 1; n <= N; n++) gel(w,n) = cgetr(prec);
    1325         266 :   av0 = avma;
    1326         266 :   if (N < Nmin) Nmin = N;
    1327         266 :   if (!eC) eC = mpexp(C);
    1328         266 :   en = eC; affrr(eint1p(C, en), gel(w,1));
    1329        3920 :   for (n = 2; n <= Nmin; n++)
    1330             :   {
    1331             :     pari_sp av2;
    1332        3654 :     en = mulrr(en,eC); /* exp(n C) */
    1333        3654 :     av2 = avma;
    1334        3654 :     affrr(eint1p(mulru(C,n), en), gel(w,n));
    1335        3654 :     set_avma(av2);
    1336             :   }
    1337         266 :   if (Nmin == N) { set_avma(av0); return w; }
    1338             : 
    1339         259 :   DL = prec2nbits_mul(prec, M_LN2) + 5;
    1340         259 :   jmin = ceil(DL/log((double)N)) + 1;
    1341         259 :   jmax = ceil(DL/log((double)Nmin)) + 1;
    1342         259 :   v = sum_jall(C, jmax, prec);
    1343         259 :   en = powrs(eC, -N); /* exp(-N C) */
    1344         259 :   affrr(eint1p(mulru(C,N), invr(en)), gel(w,N));
    1345        7038 :   for (j = jmin, n = N-1; j <= jmax; j++)
    1346             :   {
    1347        6779 :     long limN = maxss((long)ceil(exp(DL/j)), Nmin);
    1348             :     GEN polsh;
    1349        6779 :     setlg(v, j+1);
    1350        6779 :     polsh = RgV_to_RgX_reverse(v, 0);
    1351      488610 :     for (; n >= limN; n--)
    1352             :     {
    1353      481831 :       pari_sp av2 = avma;
    1354      481831 :       GEN S = divri(mulrr(en, rX_s_eval(polsh, -n)), powuu(n,j));
    1355             :       /* w[n+1] - exp(-n C) * polsh(-n) / (-n)^j */
    1356      481831 :       GEN c = odd(j)? addrr(gel(w,n+1), S) : subrr(gel(w,n+1), S);
    1357      481831 :       affrr(c, gel(w,n)); set_avma(av2);
    1358      481831 :       en = mulrr(en,eC); /* exp(-n C) */
    1359             :     }
    1360             :   }
    1361         259 :   set_avma(av0); return w;
    1362             : }
    1363             : 
    1364             : /* erfc via numerical integration : assume real(x)>=1 */
    1365             : static GEN
    1366          14 : cxerfc_r1(GEN x, long prec)
    1367             : {
    1368             :   GEN h, h2, eh2, denom, res, lambda;
    1369             :   long u, v;
    1370          14 :   const double D = prec2nbits_mul(prec, M_LN2);
    1371          14 :   const long npoints = (long)ceil(D/M_PI)+1;
    1372          14 :   pari_sp av = avma;
    1373             :   {
    1374          14 :     double t = exp(-2*M_PI*M_PI/D); /* ~exp(-2*h^2) */
    1375          14 :     v = 30; /* bits that fit in both long and double mantissa */
    1376          14 :     u = (long)floor(t*(1L<<v));
    1377             :     /* define exp(-2*h^2) to be u*2^(-v) */
    1378             :   }
    1379          14 :   incrprec(prec);
    1380          14 :   x = gtofp(x,prec);
    1381          14 :   eh2 = sqrtr_abs(rtor(shiftr(dbltor(u),-v),prec));
    1382          14 :   h2 = negr(logr_abs(eh2));
    1383          14 :   h = sqrtr_abs(h2);
    1384          14 :   lambda = gdiv(x,h);
    1385          14 :   denom = gsqr(lambda);
    1386             :   { /* res = h/x + 2*x*h*sum(k=1,npoints,exp(-(k*h)^2)/(lambda^2+k^2)); */
    1387             :     GEN Uk; /* = exp(-(kh)^2) */
    1388          14 :     GEN Vk = eh2;/* = exp(-(2k+1)h^2) */
    1389          14 :     pari_sp av2 = avma;
    1390             :     long k;
    1391             :     /* k = 0 moved out for efficiency */
    1392          14 :     denom = gaddsg(1,denom);
    1393          14 :     Uk = Vk;
    1394          14 :     Vk = mulur(u,Vk); shiftr_inplace(Vk, -v);
    1395          14 :     res = gdiv(Uk, denom);
    1396         420 :     for (k = 1; k < npoints; k++)
    1397             :     {
    1398         406 :       if ((k & 255) == 0) gerepileall(av2,4,&denom,&Uk,&Vk,&res);
    1399         406 :       denom = gaddsg(2*k+1,denom);
    1400         406 :       Uk = mpmul(Uk,Vk);
    1401         406 :       Vk = mulur(u,Vk); shiftr_inplace(Vk, -v);
    1402         406 :       res = gadd(res, gdiv(Uk, denom));
    1403             :     }
    1404             :   }
    1405          14 :   res = gmul(res, gshift(lambda,1));
    1406             :   /* 0 term : */
    1407          14 :   res = gadd(res, ginv(lambda));
    1408          14 :   res = gmul(res, gdiv(gexp(gneg(gsqr(x)), prec), mppi(prec)));
    1409          14 :   if (rtodbl(real_i(x)) < sqrt(D))
    1410             :   {
    1411          14 :     GEN t = gmul(divrr(Pi2n(1,prec),h), x);
    1412          14 :     res = gsub(res, gdivsg(2, cxexpm1(t, prec)));
    1413             :   }
    1414          14 :   return gerepileupto(av,res);
    1415             : }
    1416             : 
    1417             : GEN
    1418          35 : gerfc(GEN x, long prec)
    1419             : {
    1420             :   GEN z, xr, xi, res;
    1421             :   pari_sp av;
    1422             : 
    1423          35 :   x = trans_fix_arg(&prec,&x,&xr,&xi, &av,&res);
    1424          35 :   if (signe(xr) >= 0) {
    1425          28 :     if (cmprs(xr, 1) > 0) /* use numerical integration */
    1426          14 :       z = cxerfc_r1(x, prec);
    1427             :     else
    1428             :     { /* erfc(x) = incgam(1/2,x^2)/sqrt(Pi) */
    1429          14 :       GEN sqrtpi = sqrtr(mppi(prec));
    1430          14 :       z = incgam0(ghalf, gsqr(x), sqrtpi, prec);
    1431          14 :       z = gdiv(z, sqrtpi);
    1432             :     }
    1433             :   }
    1434             :   else
    1435             :   { /* erfc(-x)=2-erfc(x) */
    1436             :     /* FIXME could decrease prec
    1437             :     long size = nbits2extraprec((imag(x)^2-real(x)^2)/log(2));
    1438             :     prec = size > 0 ? prec : prec + size;
    1439             :     */
    1440             :     /* NOT gsubsg(2, ...) : would create a result of
    1441             :      * huge accuracy if re(x)>>1, rounded to 2 by subsequent affc_fixlg... */
    1442           7 :     z = gsub(real2n(1,prec+EXTRAPREC), gerfc(gneg(x), prec));
    1443             :   }
    1444          35 :   set_avma(av); return affc_fixlg(z, res);
    1445             : }
    1446             : 
    1447             : /***********************************************************************/
    1448             : /**                                                                   **/
    1449             : /**                      RIEMANN ZETA FUNCTION                        **/
    1450             : /**                                                                   **/
    1451             : /***********************************************************************/
    1452             : static const double log2PI = 1.83787706641;
    1453             : 
    1454             : static double
    1455        4361 : get_xinf(double beta)
    1456             : {
    1457        4361 :   const double maxbeta = 0.06415003; /* 3^(-2.5) */
    1458             :   double x0, y0, x1;
    1459             : 
    1460        4361 :   if (beta < maxbeta) return beta + pow(3*beta, 1.0/3.0);
    1461        4361 :   x0 = beta + M_PI/2.0;
    1462             :   for(;;)
    1463             :   {
    1464        9849 :     y0 = x0*x0;
    1465        7105 :     x1 = (beta+atan(x0)) * (1+y0) / y0 - 1/x0;
    1466        7105 :     if (0.99*x0 < x1) return x1;
    1467        2744 :     x0 = x1;
    1468             :   }
    1469             : }
    1470             : /* optimize for zeta( s + it, prec ), assume |s-1| > 0.1
    1471             :  * (if gexpo(u = s-1) < -5, we use the functional equation s->1-s) */
    1472             : static void
    1473        8358 : optim_zeta(GEN S, long prec, long *pp, long *pn)
    1474             : {
    1475             :   double s, t, alpha, beta, n, B;
    1476             :   long p;
    1477        8358 :   if (typ(S) == t_REAL) {
    1478         266 :     s = rtodbl(S);
    1479         266 :     t = 0.;
    1480             :   } else {
    1481        8092 :     s = rtodbl(gel(S,1));
    1482        8092 :     t = fabs( rtodbl(gel(S,2)) );
    1483             :   }
    1484             : 
    1485        8358 :   B = prec2nbits_mul(prec, M_LN2);
    1486        8358 :   if (s > 0 && !t) /* positive real input */
    1487             :   {
    1488         245 :     beta = B + 0.61 + s*(log2PI - log(s));
    1489         490 :     if (beta > 0)
    1490             :     {
    1491         238 :       p = (long)ceil(beta / 2.0);
    1492         238 :       n = fabs(s + 2*p-1)/(2*M_PI);
    1493             :     }
    1494             :     else
    1495             :     {
    1496           7 :       p = 0;
    1497           7 :       n = exp((B - M_LN2) / s);
    1498             :     }
    1499             :   }
    1500        8113 :   else if (s <= 0 || t < 0.01) /* s < 0 may occur if s ~ 0 */
    1501        3745 :   { /* TODO: the crude bounds below are generally valid. Optimize ? */
    1502        3745 :     double l,l2, la = 1.; /* heuristic */
    1503        3745 :     double rlog, ilog; dcxlog(s-1,t, &rlog,&ilog);
    1504        3745 :     l2 = (s - 0.5)*rlog - t*ilog; /* = Re( (S - 1/2) log (S-1) ) */
    1505        3745 :     l = (B - l2 + s*log2PI) / (2. * (1.+ log((double)la)));
    1506        3745 :     l2 = dabs(s, t)/2;
    1507        3745 :     if (l < l2) l = l2;
    1508        3745 :     p = (long) ceil(l); if (p < 2) p = 2;
    1509        3745 :     n = 1 + dabs(p+s/2.-.25, t/2) * la / M_PI;
    1510             :   }
    1511             :   else
    1512             :   {
    1513        4368 :     double sn = dabs(s, t), L = log(sn/s);
    1514        4368 :     alpha = B - 0.39 + L + s*(log2PI - log(sn));
    1515        4368 :     beta = (alpha+s)/t - atan(s/t);
    1516        4368 :     p = 0;
    1517        4368 :     if (beta > 0)
    1518             :     {
    1519        4361 :       beta = 1.0 - s + t * get_xinf(beta);
    1520        4361 :       if (beta > 0) p = (long)ceil(beta / 2.0);
    1521             :     }
    1522             :     else
    1523           7 :       if (s < 1.0) p = 1;
    1524        4368 :     n = p? dabs(s + 2*p-1, t) / (2*M_PI) : exp((B-M_LN2+L) / s);
    1525             :   }
    1526        8358 :   *pp = p;
    1527        8358 :   *pn = (long)ceil(n);
    1528        8358 :   if (*pp < 0 || *pn < 0) pari_err_OVERFLOW("zeta");
    1529        8358 : }
    1530             : 
    1531             : /* zeta(a*j+b), j=0..N-1, b>1, using sumalt. Johansonn's thesis, Algo 4.7.1 */
    1532             : static GEN
    1533         366 : veczetas(long a, long b, long N, long prec)
    1534             : {
    1535         366 :   const long n = ceil(2 + prec2nbits_mul(prec, M_LN2/1.7627));
    1536         366 :   pari_sp av = avma;
    1537         366 :   GEN c, d, z = zerovec(N);
    1538             :   long j, k;
    1539         366 :   c = d = int2n(2*n-1);
    1540       52313 :   for (k = n; k > 1; k--)
    1541             :   {
    1542       51947 :     GEN u, t = divii(d, powuu(k,b));
    1543       51947 :     if (!odd(k)) t = negi(t);
    1544       51947 :     gel(z,1) = addii(gel(z,1), t);
    1545       51947 :     u = powuu(k,a);
    1546     2059567 :     for (j = 1; j < N; j++)
    1547             :     {
    1548     2017724 :       t = divii(t,u); if (!signe(t)) break;
    1549     2007620 :       gel(z,j+1) = addii(gel(z,j+1), t);
    1550             :     }
    1551       51947 :     c = muluui(k,2*k-1,c);
    1552       51947 :     c = diviuuexact(c, 2*(n-k+1),n+k-1);
    1553       51947 :     d = addii(d,c);
    1554       51947 :     if (gc_needed(av,3))
    1555             :     {
    1556           9 :       if(DEBUGMEM>1) pari_warn(warnmem,"zetaBorwein, k = %ld", k);
    1557           9 :       gerepileall(av, 3, &c,&d,&z);
    1558             :     }
    1559             :   }
    1560             :   /* k = 1 */
    1561       13687 :   for (j = 1; j <= N; j++) gel(z,j) = addii(gel(z,j), d);
    1562         366 :   d = addiu(d, 1);
    1563       13687 :   for (j = 0, k = b - 1; j < N; j++, k += a)
    1564       13321 :     gel(z,j+1) = rdivii(shifti(gel(z,j+1), k), subii(shifti(d,k), d), prec);
    1565         366 :   return z;
    1566             : }
    1567             : /* zeta(a*j+b), j=0..N-1, b>1, using sumalt */
    1568             : GEN
    1569         224 : veczeta(GEN a, GEN b, long N, long prec)
    1570             : {
    1571         224 :   pari_sp av = avma;
    1572             :   long n, j, k;
    1573         224 :   GEN L, c, d, z = zerovec(N);
    1574         224 :   if (typ(a) == t_INT && typ(b) == t_INT)
    1575         126 :     return gerepilecopy(av, veczetas(itos(a),  itos(b), N, prec));
    1576          98 :   n = ceil(2 + prec2nbits_mul(prec, M_LN2/1.7627));
    1577          98 :   c = d = int2n(2*n-1);
    1578       16072 :   for (k = n; k; k--)
    1579             :   {
    1580             :     GEN u, t;
    1581       15974 :     L = logr_abs(utor(k, prec)); /* log(k) */
    1582       15974 :     t = gdiv(d, gexp(gmul(b, L), prec)); /* d / k^b */
    1583       15974 :     if (!odd(k)) t = gneg(t);
    1584       15974 :     gel(z,1) = gadd(gel(z,1), t);
    1585       15974 :     u = gexp(gmul(a, L), prec);
    1586     1021885 :     for (j = 1; j < N; j++)
    1587             :     {
    1588     1011858 :       t = gdiv(t,u); if (gexpo(t) < 0) break;
    1589     1005911 :       gel(z,j+1) = gadd(gel(z,j+1), t);
    1590             :     }
    1591       15974 :     c = muluui(k,2*k-1,c);
    1592       15974 :     c = diviuuexact(c, 2*(n-k+1),n+k-1);
    1593       15974 :     d = addii(d,c);
    1594       15974 :     if (gc_needed(av,3))
    1595             :     {
    1596           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"veczeta, k = %ld", k);
    1597           0 :       gerepileall(av, 3, &c,&d,&z);
    1598             :     }
    1599             :   }
    1600          98 :   L = mplog2(prec);
    1601        6181 :   for (j = 0; j < N; j++)
    1602             :   {
    1603        6083 :     GEN u = gsubgs(gadd(b, gmulgs(a,j)), 1);
    1604        6083 :     GEN w = gexp(gmul(u, L), prec); /* 2^u */
    1605        6083 :     gel(z,j+1) = gdiv(gmul(gel(z,j+1), w), gmul(d,gsubgs(w,1)));
    1606             :   }
    1607          98 :   return gerepilecopy(av, z);
    1608             : }
    1609             : 
    1610             : GEN
    1611       10504 : constzeta(long n, long prec)
    1612             : {
    1613       10504 :   GEN o = zetazone, z;
    1614       10504 :   long l = o? lg(o): 0;
    1615             :   pari_sp av;
    1616       10504 :   if (l > n)
    1617             :   {
    1618       10378 :     long p = realprec(gel(o,1));
    1619       10378 :     if (p >= prec) return o;
    1620             :   }
    1621         240 :   n = maxss(n, l + 15);
    1622         240 :   av = avma; z = veczetas(1, 2, n-1, prec);
    1623         240 :   zetazone = gclone(vec_prepend(z, mpeuler(prec)));
    1624         240 :   set_avma(av); guncloneNULL(o); return zetazone;
    1625             : }
    1626             : 
    1627             : /* zeta(s) using sumalt, case h=0,N=1. Assume s > 1 */
    1628             : static GEN
    1629         584 : zetaBorwein(long s, long prec)
    1630             : {
    1631         584 :   pari_sp av = avma;
    1632         584 :   const long n = ceil(2 + prec2nbits_mul(prec, M_LN2/1.7627));
    1633             :   long k;
    1634         584 :   GEN c, d, z = gen_0;
    1635         584 :   c = d = int2n(2*n-1);
    1636      134741 :   for (k = n; k; k--)
    1637             :   {
    1638      134157 :     GEN t = divii(d, powuu(k,s));
    1639      134157 :     z = odd(k)? addii(z,t): subii(z,t);
    1640      134157 :     c = muluui(k,2*k-1,c);
    1641      134157 :     c = diviuuexact(c, 2*(n-k+1),n+k-1);
    1642      134157 :     d = addii(d,c);
    1643      134157 :     if (gc_needed(av,3))
    1644             :     {
    1645           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"zetaBorwein, k = %ld", k);
    1646           0 :       gerepileall(av, 3, &c,&d,&z);
    1647             :     }
    1648             :   }
    1649         584 :   return rdivii(shifti(z, s-1), subii(shifti(d,s-1), d), prec);
    1650             : }
    1651             : 
    1652             : /* assume k != 1 */
    1653             : GEN
    1654        4718 : szeta(long k, long prec)
    1655             : {
    1656        4718 :   pari_sp av = avma;
    1657             :   GEN z;
    1658             : 
    1659        4718 :   if (!k) { z = real2n(-1, prec); setsigne(z,-1); return z; }
    1660        4711 :   if (k < 0)
    1661             :   {
    1662         252 :     if ((k&1) == 0) return gen_0;
    1663             :     /* the one value such that k < 0 and 1 - k < 0, due to overflow */
    1664         252 :     if ((ulong)k == (HIGHBIT | 1))
    1665           0 :       pari_err_OVERFLOW("zeta [large negative argument]");
    1666         252 :     k = 1-k;
    1667         252 :     z = bernreal(k, prec); togglesign(z);
    1668         252 :     return gerepileuptoleaf(av, divru(z, k));
    1669             :   }
    1670             :   /* k > 1 */
    1671        4459 :   if (k > prec2nbits(prec)+1) return real_1(prec);
    1672        4452 :   if (zetazone && realprec(gel(zetazone,1)) >= prec && lg(zetazone) > k)
    1673         140 :     return rtor(gel(zetazone, k), prec);
    1674        4312 :   if (!odd(k))
    1675             :   {
    1676             :     GEN B;
    1677        2555 :     if (!bernzone) constbern(0);
    1678        2555 :     if (k < lg(bernzone))
    1679        2100 :       B = gel(bernzone, k>>1);
    1680             :     else
    1681             :     {
    1682         455 :       if (bernbitprec(k) > prec2nbits(prec))
    1683           0 :         return gerepileupto(av, invr(inv_szeta_euler(k, prec)));
    1684         455 :       B = bernfrac(k);
    1685             :     }
    1686             :     /* B = B_k */
    1687        2555 :     z = gmul(powru(Pi2n(1, prec + EXTRAPRECWORD), k), B);
    1688        2555 :     z = divrr(z, mpfactr(k,prec));
    1689        2555 :     setsigne(z, 1); shiftr_inplace(z, -1);
    1690             :   }
    1691             :   else
    1692             :   {
    1693        1757 :     double p = prec2nbits_mul(prec,0.393); /* bit / log_2(3+sqrt(8)) */
    1694        1757 :     p = log2(p * log(p));
    1695        3514 :     z = (p * k > prec2nbits(prec))? invr(inv_szeta_euler(k, prec))
    1696        1757 :                                   : zetaBorwein(k, prec);
    1697             :   }
    1698        4312 :   return gerepileuptoleaf(av, z);
    1699             : }
    1700             : 
    1701             : /* s0 a t_INT, t_REAL or t_COMPLEX.
    1702             :  * If a t_INT, assume it's not a trivial case (i.e we have s0 > 1, odd)
    1703             :  * */
    1704             : static GEN
    1705        8379 : czeta(GEN s0, long prec)
    1706             : {
    1707        8379 :   GEN s, u, y, res, tes, sig, tau, invn2, ns, Ns, funeq_factor = NULL;
    1708             :   long i, nn, lim, lim2;
    1709        8379 :   pari_sp av0 = avma, av, av2;
    1710             :   pari_timer T;
    1711             : 
    1712        8379 :   if (DEBUGLEVEL>2) timer_start(&T);
    1713        8379 :   s = trans_fix_arg(&prec,&s0,&sig,&tau,&av,&res);
    1714        8379 :   if (typ(s0) == t_INT) return gerepileupto(av0, gzeta(s0, prec));
    1715        8372 :   if (!signe(tau)) /* real */
    1716             :   {
    1717         273 :     long e = expo(sig);
    1718         273 :     if (e >= -5 && (signe(sig) <= 0 || e < -1))
    1719             :     { /* s < 1/2 */
    1720          14 :       s = subsr(1, s);
    1721          14 :       funeq_factor = gen_1;
    1722             :     }
    1723             :   }
    1724             :   else
    1725             :   {
    1726        8099 :     u = gsubsg(1, s); /* temp */
    1727        8099 :     if (gexpo(u) < -5 || ((signe(sig) <= 0 || expo(sig) < -1) && gexpo(s) > -5))
    1728             :     { /* |1-s| < 1/32  || (real(s) < 1/2 && |imag(s)| > 1/32) */
    1729        1477 :       s = u;
    1730        1477 :       funeq_factor = gen_1;
    1731             :     }
    1732             :   }
    1733             : 
    1734        8372 :   if (funeq_factor)
    1735             :   { /* s <--> 1-s */
    1736             :     GEN t;
    1737        1491 :     sig = real_i(s);
    1738             :     /* Gamma(s) (2Pi)^-s 2 cos(Pi s/2) */
    1739        1491 :     t = gmul(ggamma(s,prec), gpow(Pi2n(1,prec), gsubgs(s0,1), prec));
    1740        1491 :     funeq_factor = gmul2n(gmul(t, gsin(gmul(Pi2n(-1,prec),s0), prec)), 1);
    1741             :   }
    1742        8372 :   if (gcmpgs(sig, prec2nbits(prec) + 1) > 0) { /* zeta(s) = 1 */
    1743          14 :     if (!funeq_factor) { set_avma(av0); return real_1(prec); }
    1744           7 :     return gerepileupto(av0, funeq_factor);
    1745             :   }
    1746        8358 :   optim_zeta(s, prec, &lim, &nn);
    1747        8358 :   if (DEBUGLEVEL>2) err_printf("lim, nn: [%ld, %ld]\n", lim, nn);
    1748        8358 :   incrprec(prec); /* one extra word of precision */
    1749             : 
    1750        8358 :   av2 = avma;
    1751        8358 :   Ns = vecpowug(nn, gneg(s), prec);
    1752        8358 :   ns = gel(Ns,nn);
    1753      263388 :   y = gmul2n(ns, -1); for (i = nn-1; i; i--) y = gadd(y, gel(Ns,i));
    1754        8358 :   if (DEBUGLEVEL>2) timer_printf(&T,"sum from 1 to N-1");
    1755        8358 :   gerepileall(av2, 2, &y,&ns);
    1756             : 
    1757        8358 :   invn2 = divri(real_1(prec), sqru(nn)); lim2 = lim<<1;
    1758        8358 :   constbern(lim);
    1759        8358 :   tes = bernfrac(lim2);
    1760             :   {
    1761             :     GEN s1, s2, s3, s4, s5;
    1762        8358 :     s2 = gmul(s, gsubgs(s,1));
    1763        8358 :     s3 = gmul2n(invn2,3);
    1764        8358 :     av2 = avma;
    1765        8358 :     s1 = gsubgs(gmul2n(s,1), 1);
    1766        8358 :     s4 = gmul(invn2, gmul2n(gaddsg(4*lim-2,s1),1));
    1767        8358 :     s5 = gmul(invn2, gadd(s2, gmulsg(lim2, gaddgs(s1, lim2))));
    1768      787441 :     for (i = lim2-2; i>=2; i -= 2)
    1769             :     {
    1770      779083 :       s5 = gsub(s5, s4);
    1771      779083 :       s4 = gsub(s4, s3);
    1772      779083 :       tes = gadd(bernfrac(i), divgunu(gmul(s5,tes), i+1));
    1773      779083 :       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        8358 :     u = gmul(gmul(tes,invn2), gmul2n(s2, -1));
    1780        8358 :     tes = gmulsg(nn, gaddsg(1, u));
    1781             :   }
    1782        8358 :   if (DEBUGLEVEL>2) timer_printf(&T,"Bernoulli sum");
    1783             :   /* y += tes n^(-s) / (s-1) */
    1784        8358 :   y = gadd(y, gmul(tes, gdiv(ns, gsubgs(s,1))));
    1785        8358 :   if (funeq_factor) y = gmul(y, funeq_factor);
    1786        8358 :   set_avma(av); return affc_fixlg(y,res);
    1787             : }
    1788             : 
    1789             : #if 0
    1790             : /* return P mod x^n where P is polynomial in x */
    1791             : static GEN
    1792             : pol_mod_xn(GEN P, long n)
    1793             : {
    1794             :   long j, l = lg(P), N = n+2;
    1795             :   GEN R;
    1796             :   if (l > N) l = N;
    1797             :   R = cgetg(N, t_POL); R[1] = evalvarn(0);
    1798             :   for (j = 2; j < l; j++) gel(R,j) = gel(P,j);
    1799             :   return normalizepol_lg(R, n+2);
    1800             : }
    1801             : 
    1802             : /* compute the values of the twisted partial
    1803             :    zeta function Z_f(a, c, s) for a in va */
    1804             : GEN
    1805             : twistpartialzeta(GEN q, long f, long c, GEN va, GEN cff)
    1806             : {
    1807             :   long j, k, lva = lg(va)-1, N = lg(cff)-1;
    1808             :   pari_sp av, av2;
    1809             :   GEN Ax, Cx, Bx, Dx, x = pol_x(0), y = pol_x(1);
    1810             :   GEN cyc, psm, rep, eta, etaf;
    1811             : 
    1812             :   cyc = gdiv(gsubgs(gpowgs(y, c), 1), gsubgs(y, 1));
    1813             :   psm = polsym(cyc, degpol(cyc) - 1);
    1814             :   eta = mkpolmod(y, cyc);
    1815             :   etaf = gpowgs(eta,f);
    1816             :   av = avma;
    1817             :   Ax  = gsubgs(gpowgs(gaddgs(x, 1), f), 1);
    1818             :   Ax  = gdiv(gmul(Ax, etaf), gsubsg(1, etaf));
    1819             :   Ax  = gerepileupto(av, RgX_to_FqX(Ax, cyc, q));
    1820             :   Cx  = Ax;
    1821             :   Bx  = gen_1;
    1822             :   av  = avma;
    1823             :   for (j = 2; j <= N; j++)
    1824             :   {
    1825             :     Bx = gadd(Bx, Cx);
    1826             :     Bx = FpXQX_red(Bx, cyc, q);
    1827             :     Cx = FpXQX_mul(Cx, Ax, cyc, q);
    1828             :     Cx = pol_mod_xn(Cx, N);
    1829             :     if (gequal0(Cx)) break;
    1830             :     if (gc_needed(av, 1))
    1831             :     {
    1832             :       if(DEBUGMEM>1) pari_warn(warnmem, "twistpartialzeta (1), j = %ld/%ld", j, N);
    1833             :       gerepileall(av, 2, &Cx, &Bx);
    1834             :     }
    1835             :   }
    1836             :   Bx  = lift_shallow(gmul(ginv(gsubsg(1, etaf)), Bx));
    1837             :   Bx  = gerepileupto(av, RgX_to_FqX(Bx, cyc, q));
    1838             :   Cx = lift_shallow(gmul(eta, gaddsg(1, x)));
    1839             :   Dx = pol_1(varn(x));
    1840             :   av2 = avma;
    1841             :   for (j = lva; j > 1; j--)
    1842             :   {
    1843             :     GEN Ex;
    1844             :     long e = va[j] - va[j-1];
    1845             :     if (e == 1)
    1846             :       Ex = Cx;
    1847             :     else
    1848             :       /* e is very small in general and actually very rarely different
    1849             :          to 1, it is always 1 for zetap (so it should be OK not to store
    1850             :          them or to compute them in a smart way) */
    1851             :       Ex = gpowgs(Cx, e);
    1852             :     Dx = gaddsg(1, FpXQX_mul(Dx, Ex, cyc, q));
    1853             :     if (gc_needed(av2, 1))
    1854             :     {
    1855             :       if(DEBUGMEM>1)
    1856             :         pari_warn(warnmem, "twistpartialzeta (2), j = %ld/%ld", lva-j, lva);
    1857             :       Dx = gerepileupto(av2, FpXQX_red(Dx, cyc, q));
    1858             :     }
    1859             :   }
    1860             :   Dx = FpXQX_mul(Dx, Cx, cyc, q); /* va[1] = 1 */
    1861             :   Bx = gerepileupto(av, FpXQX_mul(Dx, Bx, cyc, q));
    1862             :   rep = gen_0;
    1863             :   av2 = avma;
    1864             :   for (k = 1; k <= N; k++)
    1865             :   {
    1866             :     GEN p2, ak = polcoef_i(Bx, k, 0);
    1867             :     p2  = quicktrace(ak, psm);
    1868             :     rep = modii(addii(rep, mulii(gel(cff, k), p2)), q);
    1869             :     if (gc_needed(av2, 1))
    1870             :     {
    1871             :       if(DEBUGMEM>1) pari_warn(warnmem, "twistpartialzeta (3), j = %ld/%ld", k, N);
    1872             :       rep = gerepileupto(av2, rep);
    1873             :     }
    1874             :   }
    1875             :   return rep;
    1876             : }
    1877             : 
    1878             : /* initialize the roots of unity for the computation
    1879             :    of the Teichmuller character (also the values of f and c) */
    1880             : GEN
    1881             : init_teich(ulong p, GEN q, long prec)
    1882             : {
    1883             :   GEN vz, gp = utoipos(p);
    1884             :   pari_sp av = avma;
    1885             :   long j;
    1886             : 
    1887             :   if (p == 2UL)
    1888             :     return NULL;
    1889             :   else
    1890             :   { /* primitive (p-1)-th root of 1 */
    1891             :     GEN z, z0 = Zp_sqrtnlift(gen_1, utoipos(p-1), pgener_Fp(gp), gp, prec);
    1892             :     z = z0;
    1893             :     vz = cgetg(p, t_VEC);
    1894             :     for (j = 1; j < (long)p-2; j++)
    1895             :     {
    1896             :       gel(vz, umodiu(z, p)) = z; /* z = z0^i */
    1897             :       z = modii(mulii(z, z0), q);
    1898             :     }
    1899             :     gel(vz, umodiu(z, p)) = z; /* z = z0^(p-2) */
    1900             :     gel(vz,1) = gen_1; /* z0^(p-1) */
    1901             :   }
    1902             :   return gerepileupto(av, gcopy(vz));
    1903             : }
    1904             : 
    1905             : /* compute phi^(m)_s(x); s must be an integer */
    1906             : GEN
    1907             : phi_ms(ulong p, GEN q, long m, GEN s, long x, GEN vz)
    1908             : {
    1909             :   long xp = x % p;
    1910             :   GEN p1, p2;
    1911             : 
    1912             :   if (!xp) return gen_0;
    1913             :   if (vz)
    1914             :     p1 =gel(vz,xp); /* vz[x] = Teichmuller(x) */
    1915             :   else
    1916             :     p1 = (x & 2)? gen_m1: gen_1;
    1917             :   p1 = Fp_pow(p1, addis(s, m), q);
    1918             :   p2 = Fp_pow(stoi(x), negi(s), q);
    1919             :   return modii(mulii(p1,p2), q);
    1920             : }
    1921             : 
    1922             : /* compute the first N coefficients of the Mahler expansion
    1923             :    of phi^m_s skipping the first one (which is zero) */
    1924             : GEN
    1925             : coeff_of_phi_ms(ulong p, GEN q, long m, GEN s, long N, GEN vz)
    1926             : {
    1927             :   GEN qs2 = shifti(q, -1), cff = zerovec(N);
    1928             :   pari_sp av;
    1929             :   long k, j;
    1930             : 
    1931             :   av = avma;
    1932             :   for (k = 1; k <= N; k++)
    1933             :   {
    1934             :     gel(cff, k) = phi_ms(p, q, m, s, k, vz);
    1935             :     if (gc_needed(av, 2))
    1936             :     {
    1937             :       if(DEBUGMEM>1)
    1938             :         pari_warn(warnmem, "coeff_of_phi_ms (1), k = %ld/%ld", N-k, N);
    1939             :       cff = gerepileupto(av, gcopy(cff));
    1940             :     }
    1941             :   }
    1942             :   for (j = N; j > 1; j--)
    1943             :   {
    1944             :     GEN b = subii(gel(cff, j), gel(cff, j-1));
    1945             :     gel(cff, j) = centermodii(b, q, qs2);
    1946             :     if (gc_needed(av, 2))
    1947             :     {
    1948             :       if(DEBUGMEM>1)
    1949             :         pari_warn(warnmem, "coeff_of_phi_ms (2), j = %ld/%ld", N-j, N);
    1950             :       cff = gerepileupto(av, gcopy(cff));
    1951             :     }
    1952             :   }
    1953             :   for (k = 1; k < N; k++)
    1954             :     for (j = N; j > k; j--)
    1955             :     {
    1956             :       GEN b = subii(gel(cff, j), gel(cff, j-1));
    1957             :       gel(cff, j) = centermodii(b, q, qs2);
    1958             :       if (gc_needed(av, 2))
    1959             :       {
    1960             :         if(DEBUGMEM>1)
    1961             :           pari_warn(warnmem, "coeff_of_phi_ms (3), (k,j) = (%ld,%ld)/%ld",
    1962             :               k, N-j, N);
    1963             :         cff = gerepileupto(av, gcopy(cff));
    1964             :       }
    1965             :     }
    1966             :   k = N; while(gequal0(gel(cff, k))) k--;
    1967             :   setlg(cff, k+1);
    1968             :   if (DEBUGLEVEL > 2)
    1969             :     err_printf("  coeff_of_phi_ms: %ld coefficients kept out of %ld\n",
    1970             :                k, N);
    1971             :   return gerepileupto(av, cff);
    1972             : }
    1973             : 
    1974             : static long
    1975             : valfact(long N, ulong p)
    1976             : {
    1977             :   long f = 0;
    1978             :   while (N > 1)
    1979             :   {
    1980             :     N /= p;
    1981             :     f += N;
    1982             :   }
    1983             :   return f;
    1984             : }
    1985             : 
    1986             : static long
    1987             : number_of_terms(ulong p, long prec)
    1988             : {
    1989             :   long N, f;
    1990             : 
    1991             :   if (prec == 0) return p;
    1992             :   N = (long)((p-1)*prec + (p>>1)*(log2(prec)/log2(p)));
    1993             :   N = p*(N/p);
    1994             :   f = valfact(N, p);
    1995             :   while (f > prec)
    1996             :   {
    1997             :     N = p*(N/p) - 1;
    1998             :     f -= u_lval(N+1, p);
    1999             :   }
    2000             :   while (f < prec)
    2001             :   {
    2002             :     N = p*(N/p+1);
    2003             :     f += u_lval(N, p);
    2004             :   }
    2005             :   return N;
    2006             : }
    2007             : 
    2008             : static GEN
    2009             : zetap(GEN s)
    2010             : {
    2011             :   ulong p;
    2012             :   long N, f, c, prec = precp(s);
    2013             :   pari_sp av = avma;
    2014             :   GEN gp, q, vz, is, cff, val, va, cft;
    2015             : 
    2016             :   if (valp(s) < 0) pari_err_DOMAIN("zetap", "v_p(s)", "<", gen_0, s);
    2017             :   if (!prec) prec = 1;
    2018             : 
    2019             :   gp = gel(s,2); p = itou(gp);
    2020             :   is = gtrunc(s);  /* make s an integer */
    2021             : 
    2022             :   N  = number_of_terms(p, prec);
    2023             :   q  = powiu(gp, prec);
    2024             : 
    2025             :   /* initialize the roots of unity for the computation
    2026             :      of the Teichmuller character (also the values of f and c) */
    2027             :   if (DEBUGLEVEL > 1) err_printf("zetap: computing (p-1)th roots of 1\n");
    2028             :   vz = init_teich(p, q, prec);
    2029             :   if (p == 2UL) {  f = 4; c = 3; } else { f = (long)p; c = 2; }
    2030             : 
    2031             :   /* compute the first N coefficients of the Mahler expansion
    2032             :      of phi^(-1)_s skipping the first one (which is zero) */
    2033             :   if (DEBUGLEVEL > 1)
    2034             :     err_printf("zetap: computing Mahler expansion of phi^(-1)_s\n");
    2035             :   cff = coeff_of_phi_ms(p, q, -1, is, N, vz);
    2036             : 
    2037             :   /* compute the coefficients of the power series corresponding
    2038             :      to the twisted partial zeta function Z_f(a, c, s) for a in va */
    2039             :   /* The line below looks a bit stupid but it is to keep the
    2040             :      possibility of later adding p-adic Dirichlet L-functions */
    2041             :   va = identity_perm(f - 1);
    2042             :   if (DEBUGLEVEL > 1)
    2043             :     err_printf("zetap: computing values of twisted partial zeta functions\n");
    2044             :   val = twistpartialzeta(q, f, c, va, cff);
    2045             : 
    2046             :   /* sum over all a's the coefficients of the twisted
    2047             :      partial zeta functions and integrate */
    2048             :   if (DEBUGLEVEL > 1)
    2049             :     err_printf("zetap: multiplying by correcting factor\n");
    2050             : 
    2051             :   /* multiply by the corrective factor */
    2052             :   cft = gsubgs(gmulsg(c, phi_ms(p, q, -1, is, c, vz)), 1);
    2053             :   val = gdiv(val, cft);
    2054             : 
    2055             :   /* adjust the precision and return */
    2056             :   return gerepileupto(av, cvtop(val, gp, prec));
    2057             : }
    2058             : #else
    2059             : /* s1 = s-1 or NULL (if s=1) */
    2060             : static GEN
    2061         168 : hurwitzp_i(GEN cache, GEN s1, GEN x, GEN p, long prec)
    2062             : {
    2063         168 :   long j, J = lg(cache) - 2;
    2064             :   GEN S, x2, x2j;
    2065             : 
    2066         168 :   x = ginv(gadd(x, zeropadic_shallow(p, prec)));
    2067         168 :   S = s1? gmul2n(gmul(s1, x), -1): gadd(Qp_log(x), gmul2n(x, -1));
    2068         168 :   x2j = x2 = gsqr(x); S = gaddgs(S,1);
    2069         168 :   for (j = 1;; j++)
    2070             :   {
    2071         511 :     S = gadd(S, gmul(gel(cache, j+1), x2j));
    2072         511 :     if (j == J) break;
    2073         343 :     x2j = gmul(x2, x2j);
    2074             :   }
    2075         168 :   if (s1) S = gmul(gdiv(S, s1), Qp_exp(gmul(s1, Qp_log(x))));
    2076         168 :   return S;
    2077             : }
    2078             : 
    2079             : static GEN
    2080         161 : init_cache(long prec, long p, GEN s)
    2081             : {
    2082         161 :   long j, fls = !gequal1(s), J = (((p==2)? (prec >> 1): prec) + 2) >> 1;
    2083         161 :   GEN C = gen_1, cache = bernvec(J);
    2084         630 :   for (j = 1; j <= J; j++)
    2085             :   { /* B_{2j} * binomial(1-s, 2j) */
    2086         469 :     GEN t = (j > 1 || fls) ? gmul(gaddgs(s, 2*j-3), gaddgs(s, 2*j-2)) : s;
    2087         469 :     C = gdiv(gmul(C, t), mulss(2*j, 2*j-1));
    2088         469 :     gel(cache, j+1) = gmul(gel(cache, j+1), C);
    2089             :   }
    2090         161 :   return cache;
    2091             : }
    2092             : 
    2093             : static GEN
    2094          14 : zetap_i(GEN s, long D)
    2095             : {
    2096          14 :   pari_sp av = avma;
    2097          14 :   GEN cache, S, s1, gm, va, gp = gel(s,2);
    2098          14 :   ulong a, p = itou(gp), m;
    2099          14 :   long prec = valp(s) + precp(s);
    2100             : 
    2101          14 :   if (D <= 0) pari_err_DOMAIN("p-adic L-function", "D", "<=", gen_0, stoi(D));
    2102          14 :   if (!uposisfundamental(D))
    2103           0 :     pari_err_TYPE("p-adic L-function [D not fundamental]", stoi(D));
    2104          14 :   if (D == 1 && gequal1(s))
    2105           0 :     pari_err_DOMAIN("p-adic zeta", "argument", "=", gen_1, s);
    2106          14 :   if (prec <= 0) prec = 1;
    2107          14 :   cache = init_cache(prec, p, s);
    2108          14 :   m = ulcm(D, p == 2? 4: p);
    2109          14 :   gm = stoi(m);
    2110          14 :   va = coprimes_zv(m);
    2111          14 :   S = gen_0; s1 = gsubgs(s,1); if (gequal0(s1)) s1 = NULL;
    2112          42 :   for (a = 1; a <= (m >> 1); a++)
    2113          28 :     if (va[a])
    2114             :     {
    2115          21 :       GEN z = hurwitzp_i(cache, s1, gdivsg(a,gm), gp, prec);
    2116          21 :       S = gadd(S, gmulsg(kross(D,a), z));
    2117             :     }
    2118          14 :   S = gdiv(gmul2n(S, 1), gm);
    2119          14 :   if (D > 1)
    2120             :   {
    2121           0 :     gm = gadd(gm, zeropadic_shallow(gp, prec));
    2122           0 :     S = gmul(S, Qp_exp(gmul(gsubsg(1, s), Qp_log(gm))));
    2123             :   }
    2124          14 :   return gerepileupto(av, S);
    2125             : }
    2126             : static GEN
    2127          14 : zetap(GEN s) { return zetap_i(s, 1); }
    2128             : #endif
    2129             : 
    2130             : /* v a t_VEC/t_COL */
    2131             : int
    2132          14 : RgV_is_arithprog(GEN v, GEN *a, GEN *b)
    2133             : {
    2134          14 :   pari_sp av = avma, av2;
    2135          14 :   long i, n = lg(v)-1;
    2136          14 :   if (n == 0) { *a = *b = gen_0; return 1; }
    2137          14 :   *a = gel(v,1);
    2138          14 :   if (n == 1) { * b = gen_0; return 1; }
    2139          14 :   *b = gsub(gel(v,2), *a); av2 = avma;
    2140          28 :   for (i = 2; i < n; i++)
    2141          14 :     if (!gequal(*b, gsub(gel(v,i+1), gel(v,i)))) return gc_int(av,0);
    2142          14 :   return gc_int(av2,1);
    2143             : }
    2144             : 
    2145             : GEN
    2146       13090 : gzeta(GEN x, long prec)
    2147             : {
    2148       13090 :   pari_sp av = avma;
    2149             :   GEN y;
    2150       13090 :   if (gequal1(x)) pari_err_DOMAIN("zeta", "argument", "=", gen_1, x);
    2151       13083 :   switch(typ(x))
    2152             :   {
    2153        4487 :     case t_INT:
    2154        4487 :       if (is_bigint(x))
    2155             :       {
    2156          21 :         if (signe(x) > 0) return real_1(prec);
    2157          14 :         if (mod2(x) == 0) return real_0(prec);
    2158           7 :         pari_err_OVERFLOW("zeta [large negative argument]");
    2159             :       }
    2160        4466 :       return szeta(itos(x),prec);
    2161        8379 :     case t_REAL: case t_COMPLEX: return czeta(x,prec);
    2162          14 :     case t_PADIC: return zetap(x);
    2163          21 :     case t_VEC: case t_COL:
    2164             :     {
    2165             :       GEN a, b;
    2166          21 :       long n = lg(x) - 1;
    2167          21 :       if (n > 1 && RgV_is_arithprog(x, &a, &b))
    2168             :       {
    2169          14 :         if (!is_real_t(typ(a)) || !is_real_t(typ(b)) || gcmpgs(a, 1) <= 0)
    2170           0 :         { set_avma(av); break; }
    2171          14 :         a = veczeta(b, a, n, prec);
    2172          14 :         settyp(a, typ(x)); return a;
    2173             :       }
    2174             :     }
    2175             :     default:
    2176         189 :       if (!(y = toser_i(x))) break;
    2177          21 :       return gerepileupto(av, lfun(gen_1,y,prec2nbits(prec)));
    2178             :   }
    2179         168 :   return trans_eval("zeta",gzeta,x,prec);
    2180             : }
    2181             : 
    2182             : /********************************************************/
    2183             : /*                   Hurwitz zeta function              */
    2184             : /********************************************************/
    2185             : 
    2186             : static GEN
    2187         175 : hurwitzp(GEN s, GEN x)
    2188             : {
    2189         175 :   GEN s1, T = (typ(s) == t_PADIC)? s: x, gp = gel(T,2);
    2190         175 :   long p = itou(gp), vqp = (p==2)? 2: 1, prec = maxss(1, valp(T) + precp(T));
    2191             : 
    2192         175 :   if (s == T) x = gadd(x, zeropadic_shallow(gp, prec));
    2193           7 :   else        s = gadd(s, zeropadic_shallow(gp, prec));
    2194             :   /* now both s and x are t_PADIC */
    2195         175 :   if (valp(x) > -vqp)
    2196             :   {
    2197          28 :     GEN S = gen_0;
    2198          28 :     long j, M = (p==2)? 4: p;
    2199         189 :     for (j = 0; j < M; j++)
    2200             :     {
    2201         161 :       GEN y = gaddsg(j, x);
    2202         161 :       if (valp(y) <= 0) S = gadd(S, hurwitzp(s, gdivgs(y, M)));
    2203             :     }
    2204          28 :     return gdivgs(S, M);
    2205             :   }
    2206         147 :   if (valp(s) <= 1/(p-1) - vqp)
    2207           0 :     pari_err_DOMAIN("hurwitzp", "v(s)", "<=", stoi(1/(p-1)-vqp), s);
    2208         147 :   s1 = gsubgs(s,1); if (gequal0(s1)) s1 = NULL;
    2209         147 :   return hurwitzp_i(init_cache(prec,p,s), s1, x, gp, prec);
    2210             : }
    2211             : 
    2212             : static void
    2213        5138 : binsplit(GEN *pP, GEN *pR, GEN aN2, GEN isqaN, GEN s, long j, long k, long prec)
    2214             : {
    2215        5138 :   if (j + 1 == k)
    2216             :   {
    2217        2618 :     long j2 = j << 1;
    2218             :     GEN P;
    2219        2618 :     if (!j) P = gdiv(s, aN2);
    2220             :     else
    2221             :     {
    2222        2520 :       P = gmul(gaddgs(s, j2-1), gaddgs(s, j2));
    2223        2520 :       P = gdivgs(gmul(P, isqaN), (j2+1) * (j2+2));
    2224             :     }
    2225        2618 :     if (pP) *pP = P;
    2226        2618 :     if (pR) *pR = gmul(bernreal(j2+2, prec), P);
    2227             :   }
    2228             :   else
    2229             :   {
    2230             :     GEN P1, R1, P2, R2;
    2231        2520 :     binsplit(&P1,pR? &R1: NULL, aN2, isqaN, s, j, (j+k) >> 1, prec);
    2232        2520 :     binsplit(pP? &P2: NULL, pR? &R2: NULL, aN2, isqaN, s, (j+k) >> 1, k, prec);
    2233        2520 :     if (pP) *pP = gmul(P1,P2);
    2234        2520 :     if (pR) *pR = gadd(R1, gmul(P1, R2));
    2235             :   }
    2236        5138 : }
    2237             : 
    2238             : /* New zetahurwitz, from Fredrik Johansson. */
    2239             : GEN
    2240         231 : zetahurwitz(GEN s, GEN x, long der, long bitprec)
    2241             : {
    2242         231 :   pari_sp av = avma;
    2243         231 :   GEN al, ral, ral0, Nx, S1, S2, S3, N2, rx, sch = NULL, s0 = s, y;
    2244         231 :   long j, k, m, N, precinit = nbits2prec(bitprec), prec = precinit;
    2245         231 :   long fli = 0, v, prpr;
    2246             :   pari_timer T;
    2247             : 
    2248         231 :   if (der < 0) pari_err_DOMAIN("zetahurwitz", "der", "<", gen_0, stoi(der));
    2249         231 :   if (der)
    2250             :   {
    2251             :     GEN z;
    2252          21 :     if (!is_scalar_t(typ(s)))
    2253             :     {
    2254           7 :       z = deriv(zetahurwitz(s, x, der - 1, bitprec), -1);
    2255           7 :       z = gdiv(z, deriv(s, -1));
    2256             :     }
    2257             :     else
    2258             :     {
    2259             :       GEN sser;
    2260          14 :       if (gequal1(s)) pari_err_DOMAIN("zetahurwitz", "s", "=", gen_1, s0);
    2261          14 :       sser = gadd(gadd(s, pol_x(0)), zeroser(0, der + 2));
    2262          14 :       z = zetahurwitz(sser, x, 0, bitprec + der * log2(der));
    2263          14 :       z = gmul(mpfact(der), polcoef(z, der, -1));
    2264             :     }
    2265          21 :     return gerepileupto(av,z);
    2266             :   }
    2267         210 :   if (typ(x) == t_PADIC || typ(s) == t_PADIC)
    2268          35 :     return gerepilecopy(av, hurwitzp(s, x));
    2269         175 :   switch(typ(x))
    2270             :   {
    2271         161 :     case t_INT: case t_REAL: case t_FRAC: case t_COMPLEX:
    2272         161 :       rx = ground(real_i(x));
    2273         161 :       if (signe(rx) <= 0 && gexpo(gsub(x, rx)) < 17 - bitprec)
    2274           0 :         pari_err_DOMAIN("zetahurwitz", "x", "<=", gen_0, x);
    2275         161 :       break;
    2276          14 :     default:
    2277          14 :       if (!(y = toser_i(x))) pari_err_TYPE("zetahurwitz", x);
    2278           7 :       x = y; rx = ground(polcoef_i(x, 0, -1));
    2279           7 :       if (typ(rx) != t_INT) pari_err_TYPE("zetahurwitz", x);
    2280             :   }
    2281         168 :   switch (typ(s))
    2282             :   {
    2283          98 :     case t_INT: case t_REAL: case t_FRAC: case t_COMPLEX: break;
    2284          70 :     default:
    2285          70 :       if (!(y = toser_i(s))) pari_err_TYPE("zetahurwitz", s);
    2286          70 :       if (valp(y) < 0) pari_err_DOMAIN("zetahurwitz", "val(s)", "<", gen_0, s);
    2287          70 :       s0 = polcoef_i(y, 0, -1);
    2288          70 :       switch(typ(s0))
    2289             :       {
    2290          56 :         case t_INT: case t_REAL: case t_COMPLEX: case t_FRAC: break;
    2291           7 :         case t_PADIC: pari_err_IMPL("zetahurwitz(t_SER of t_PADIC)");
    2292           7 :         default: pari_err_TYPE("zetahurwitz", s0);
    2293             :       }
    2294          56 :       sch = gequal0(s0)? y: serchop0(y);
    2295          56 :       v = valp(sch);
    2296          56 :       prpr = (lg(y) + v + 1)/v; if (gequal1(s0)) prpr += v;
    2297          56 :       s = gadd(gadd(s0, pol_x(0)), zeroser(0, prpr));
    2298             :     }
    2299         154 :   al = gneg(s0); ral = real_i(al); ral0 = ground(ral);
    2300         154 :   if (gequal1(s0) && (!sch || gequal0(sch)))
    2301          14 :     pari_err_DOMAIN("zetahurwitz", "s", "=", gen_1, s0);
    2302         140 :   fli = (gsigne(ral0) >= 0 && gexpo(gsub(al, ral0)) < 17 - bitprec);
    2303         140 :   if (!sch && fli)
    2304             :   { /* al ~ non negative integer */
    2305           0 :     k = itos(gceil(ral)) + 1;
    2306           0 :     if (odd(k)) k++;
    2307           0 :     N = 1;
    2308             :   }
    2309             :   else
    2310             :   {
    2311         140 :     const double D = (typ(s) == t_INT)? 0.24: 0.4;
    2312         140 :     GEN C, rs = real_i(gsubsg(1, s0));
    2313         140 :     long ebit = 0;
    2314         140 :     if (fli) al = gadd(al, ghalf); /* hack */
    2315         140 :     if (gsigne(rs) > 0)
    2316             :     {
    2317          42 :       ebit = (long)(ceil(gtodouble(rs)*expu(bitprec)));
    2318          42 :       bitprec += ebit; prec = nbits2prec(bitprec);
    2319          42 :       x = gprec_w(x, prec);
    2320          42 :       s = gprec_w(s, prec);
    2321          42 :       al = gprec_w(al, prec);
    2322          42 :       if (sch) sch = gprec_w(sch, prec);
    2323             :     }
    2324         140 :     k = maxss(itos(gceil(gadd(ral, ghalf))) + 1, 50);
    2325         140 :     k = maxss(k, (long)(D*bitprec));
    2326         140 :     if (odd(k)) k++;
    2327         140 :     C = gmulsg(2, gmul(binomial(al, k+1), gdivgs(bernfrac(k+2), k+2)));
    2328         140 :     C = gmul2n(gabs(C,LOWDEFAULTPREC), bitprec);
    2329         140 :     C = gadd(gpow(C, ginv(gsubsg(k+1, ral)), LOWDEFAULTPREC),
    2330             :              gabs(gsubsg(1,rx), LOWDEFAULTPREC));
    2331         140 :     C = gceil(polcoef_i(C, 0, -1));
    2332         140 :     if (typ(C) != t_INT) pari_err_TYPE("zetahurwitz",s);
    2333         140 :     N = itos(gceil(C));
    2334         140 :     if (N < 1) N = 1;
    2335             :   }
    2336         140 :   N = maxss(N, 1 - itos(rx));
    2337         140 :   al = gneg(s);
    2338         140 :   if (DEBUGLEVEL>2) timer_start(&T);
    2339         140 :   incrprec(prec);
    2340         140 :   Nx = gmul(real_1(prec), gaddsg(N - 1, x));
    2341         140 :   S1 = S3 = gpow(Nx, al, prec);
    2342        2702 :   for (m = N - 2; m >= 0; m--) S1 = gadd(S1, gpow(gaddsg(m,x), al, prec));
    2343         140 :   if (DEBUGLEVEL>2) timer_printf(&T,"sum from 0 to N - 1");
    2344         140 :   constbern(k >> 1);
    2345         140 :   N2 = ginv(gsqr(Nx));
    2346         140 :   if (typ(s0) == t_INT)
    2347             :   {
    2348          42 :     S2 = gen_0;
    2349        1274 :     for (j = k; j >= 2; j -= 2)
    2350             :     {
    2351        1232 :       GEN t = gsubgs(al, j), u = gmul(t, gaddgs(t, 1));
    2352        1232 :       u = gmul(gdivgs(u, j*(j+1)), gmul(S2, N2));
    2353        1232 :       S2 = gadd(gdivgs(bernfrac(j), j), u);
    2354             :     }
    2355          42 :     S2 = gmul(S2, gdiv(al, Nx));
    2356             :   }
    2357             :   else
    2358             :   {
    2359          98 :     binsplit(NULL,&S2, gmul2n(Nx,1), N2, s, 0, k >> 1, prec);
    2360          98 :     S2 = gneg(S2);
    2361             :   }
    2362         140 :   S2 = gadd(ghalf, S2);
    2363         140 :   if (DEBUGLEVEL>2) timer_printf(&T,"Bernoulli sum");
    2364         140 :   S2 = gmul(S3, gadd(gdiv(Nx, gaddsg(1, al)), S2));
    2365         140 :   S1 = gprec_wtrunc(gsub(S1, S2), precinit);
    2366         140 :   if (sch) return gerepileupto(av, gsubst(S1, 0, sch));
    2367          91 :   return gerepilecopy(av, S1);
    2368             : }
    2369             : 
    2370             : /***********************************************************************/
    2371             : /**                                                                   **/
    2372             : /**                    FONCTIONS POLYLOGARITHME                       **/
    2373             : /**                                                                   **/
    2374             : /***********************************************************************/
    2375             : 
    2376             : /* returns H_n = 1 + 1/2 + ... + 1/n, as a rational number (n "small") */
    2377             : static GEN
    2378          21 : Harmonic(long n)
    2379             : {
    2380          21 :   GEN h = gen_1;
    2381             :   long i;
    2382          28 :   for (i=2; i<=n; i++) h = gadd(h, mkfrac(gen_1, utoipos(i)));
    2383          21 :   return h;
    2384             : }
    2385             : /* smallish k such that bernbitprec(K) > bit + Kdz, K = 2k+4 */
    2386             : static long
    2387          21 : get_k(double dz, long bit)
    2388             : {
    2389             :   long a, b;
    2390          21 :   for (b = 128;; b <<= 1)
    2391          21 :     if (bernbitprec(b) > bit + b*dz) break;
    2392          21 :   if (b == 128) return 128;
    2393           0 :   a = b >> 1;
    2394           0 :   while (b - a > 64)
    2395             :   {
    2396           0 :     long c = (a+b) >> 1;
    2397           0 :     if (bernbitprec(c) > bit + c*dz) b = c; else a = c;
    2398             :   }
    2399           0 :   return b >> 1;
    2400             : }
    2401             : 
    2402             : /* m >= 2. Validity domain |log x| < 2*Pi, contains log |x| < 5.44,
    2403             :  * Li_m(x = e^z) = sum_{n >= 0} zeta(m-n) z^n / n!
    2404             :  *    with zeta(1) := H_{m-1} - log(-z) */
    2405             : static GEN
    2406          21 : cxpolylog(long m, GEN x, long prec)
    2407             : {
    2408             :   long li, n, k, ksmall, real;
    2409             :   GEN vz, z, Z, h, q, s, S;
    2410             :   pari_sp av;
    2411             :   double dz;
    2412             :   pari_timer T;
    2413             : 
    2414          21 :   if (gequal1(x)) return szeta(m,prec);
    2415             :   /* x real <= 1 ==> Li_m(x) real */
    2416          21 :   real = (typ(x) == t_REAL && (expo(x) < 0 || signe(x) <= 0));
    2417             : 
    2418          21 :   vz = constzeta(m, prec);
    2419          21 :   z = glog(x,prec);
    2420             :   /* n = 0 */
    2421          21 :   q = gen_1; s = gel(vz, m);
    2422          28 :   for (n=1; n < m-1; n++)
    2423             :   {
    2424           7 :     q = gdivgs(gmul(q,z),n);
    2425           7 :     s = gadd(s, gmul(gel(vz,m-n), real? real_i(q): q));
    2426             :   }
    2427             :   /* n = m-1 */
    2428          21 :     q = gdivgs(gmul(q,z),n); /* multiply by "zeta(1)" */
    2429          21 :     h = gmul(q, gsub(Harmonic(m-1), glog(gneg_i(z),prec)));
    2430          21 :     s = gadd(s, real? real_i(h): h);
    2431             :   /* n = m */
    2432          21 :     q = gdivgs(gmul(q,z),m);
    2433          21 :     s = gadd(s, gdivgs(real? real_i(q): q, -2)); /* zeta(0) = -1/2 */
    2434             :   /* n = m+1 */
    2435          21 :     q = gdivgs(gmul(q,z),m+1); /* = z^(m+1) / (m+1)! */
    2436          21 :     s = gadd(s, gdivgs(real? real_i(q): q, -12)); /* zeta(-1) = -1/12 */
    2437             : 
    2438          21 :   li = -(prec2nbits(prec)+1);
    2439          21 :   if (DEBUGLEVEL) timer_start(&T);
    2440          21 :   dz = dbllog2(z) - log2PI; /*  ~ log2(|z|/2Pi) */
    2441             :   /* sum_{k >= 1} zeta(-1-2k) * z^(2k+m+1) / (2k+m+1)!
    2442             :    * = 2 z^(m-1) sum_{k >= 1} zeta(2k+2) * Z^(k+1) / (2k+2)..(2k+1+m)), where
    2443             :    * Z = -(z/2Pi)^2. Stop at 2k = (li - (m-1)*Lz - m) / dz, Lz = log2 |z| */
    2444             :   /* We cut the sum in two: small values of k first */
    2445          21 :   Z = gsqr(z); av = avma;
    2446          21 :   ksmall = get_k(dz, prec2nbits(prec));
    2447          21 :   constbern(ksmall);
    2448         469 :   for(k = 1; k < ksmall; k++)
    2449             :   {
    2450         469 :     GEN t = q = divgunu(gmul(q,Z), 2*k+m); /* z^(2k+m+1)/(2k+m+1)! */
    2451         469 :     if (real) t = real_i(t);
    2452         469 :     t = gmul(t, gdivgs(bernfrac(2*k+2), 2*k+2)); /* - t * zeta(1-(2k+2)) */
    2453         469 :     s = gsub(s, t);
    2454         469 :     if (gexpo(t)  < li) return s;
    2455             :     /* large values ? */
    2456         448 :     if ((k & 0x1ff) == 0) gerepileall(av, 2, &s, &q);
    2457             :   }
    2458           0 :   if (DEBUGLEVEL>2) timer_printf(&T, "polylog: small k <= %ld", k);
    2459           0 :   Z = gneg(gsqr(gdiv(z, Pi2n(1,prec))));
    2460           0 :   q = gmul(gpowgs(z, m-1), gpowgs(Z, k+1)); /* Z^(k+1) * z^(m-1) */
    2461           0 :   S = gen_0; av = avma;
    2462           0 :   for(;; k++)
    2463           0 :   {
    2464           0 :     GEN t = q;
    2465             :     long b;
    2466           0 :     if (real) t = real_i(t);
    2467           0 :     b = prec + gexpo(t) / BITS_IN_LONG; /* decrease accuracy */
    2468           0 :     if (b == 2) break;
    2469             :     /* t * zeta(2k+2) / (2k+2)..(2k+1+m) */
    2470           0 :     t = gdiv(t, mulri(inv_szeta_euler(2*k+2, b),
    2471           0 :                       mulu_interval(2*k+2, 2*k+1+m)));
    2472           0 :     S = gadd(S, t); if (gexpo(t)  < li) break;
    2473           0 :     q = gmul(q, Z);
    2474           0 :     if ((k & 0x1ff) == 0) gerepileall(av, 2, &S, &q);
    2475             :   }
    2476           0 :   if (DEBUGLEVEL>2) timer_printf(&T, "polylog: large k <= %ld", k);
    2477           0 :   return gadd(s, gmul2n(S,1));
    2478             : }
    2479             : 
    2480             : static GEN
    2481          42 : Li1(GEN x, long prec) { return gneg(glog(gsubsg(1, x), prec)); }
    2482             : 
    2483             : static GEN
    2484         203 : polylog(long m, GEN x, long prec)
    2485             : {
    2486             :   long l, e, i, G, sx;
    2487             :   pari_sp av, av1;
    2488             :   GEN X, Xn, z, p1, p2, y, res;
    2489             : 
    2490         203 :   if (m < 0) pari_err_DOMAIN("polylog", "index", "<", gen_0, stoi(m));
    2491         203 :   if (!m) return mkfrac(gen_m1,gen_2);
    2492         203 :   if (gequal0(x)) return gcopy(x);
    2493         203 :   if (m==1) { av = avma; return gerepileupto(av, Li1(x, prec)); }
    2494             : 
    2495         168 :   l = precision(x);
    2496         168 :   if (!l) l = prec; else prec = l;
    2497         168 :   res = cgetc(l); av = avma;
    2498         168 :   x = gtofp(x, l+EXTRAPRECWORD);
    2499         168 :   e = gexpo(gnorm(x));
    2500         168 :   if (!e || e == -1) {
    2501          21 :     y = cxpolylog(m,x,prec);
    2502          21 :     set_avma(av); return affc_fixlg(y, res);
    2503             :   }
    2504         147 :   X = (e > 0)? ginv(x): x;
    2505         147 :   G = -prec2nbits(l);
    2506         147 :   av1 = avma;
    2507         147 :   y = Xn = X;
    2508         147 :   for (i=2; ; i++)
    2509             :   {
    2510       68159 :     Xn = gmul(X,Xn); p2 = gdiv(Xn,powuu(i,m));
    2511       68159 :     y = gadd(y,p2);
    2512       68159 :     if (gexpo(p2) <= G) break;
    2513             : 
    2514       68012 :     if (gc_needed(av1,1))
    2515             :     {
    2516           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"polylog");
    2517           0 :       gerepileall(av1,2, &y, &Xn);
    2518             :     }
    2519             :   }
    2520         147 :   if (e < 0) { set_avma(av); return affc_fixlg(y, res); }
    2521             : 
    2522          28 :   sx = gsigne(imag_i(x));
    2523          28 :   if (!sx)
    2524             :   {
    2525          28 :     if (m&1) sx = gsigne(gsub(gen_1, real_i(x)));
    2526          21 :     else     sx = - gsigne(real_i(x));
    2527             :   }
    2528          28 :   z = divri(mppi(l), mpfact(m-1)); setsigne(z, sx);
    2529          28 :   z = mkcomplex(gen_0, z);
    2530             : 
    2531          28 :   if (m == 2)
    2532             :   { /* same formula as below, written more efficiently */
    2533          21 :     y = gneg_i(y);
    2534          21 :     if (typ(x) == t_REAL && signe(x) < 0)
    2535           7 :       p1 = logr_abs(x);
    2536             :     else
    2537          14 :       p1 = gsub(glog(x,l), z);
    2538          21 :     p1 = gmul2n(gsqr(p1), -1); /* = (log(-x))^2 / 2 */
    2539             : 
    2540          21 :     p1 = gadd(p1, divru(sqrr(mppi(l)), 6));
    2541          21 :     p1 = gneg_i(p1);
    2542             :   }
    2543             :   else
    2544             :   {
    2545           7 :     GEN logx = glog(x,l), logx2 = gsqr(logx), vz = constzeta(m, l);
    2546           7 :     p1 = mkfrac(gen_m1,gen_2);
    2547          14 :     for (i = m-2; i >= 0; i -= 2)
    2548           7 :       p1 = gadd(gel(vz, m-i), gmul(p1,gdivgs(logx2,(i+1)*(i+2))));
    2549           7 :     if (m&1) p1 = gmul(logx,p1); else y = gneg_i(y);
    2550           7 :     p1 = gadd(gmul2n(p1,1), gmul(z,gpowgs(logx,m-1)));
    2551           7 :     if (typ(x) == t_REAL && signe(x) < 0) p1 = real_i(p1);
    2552             :   }
    2553          28 :   y = gadd(y,p1);
    2554          28 :   set_avma(av); return affc_fixlg(y, res);
    2555             : }
    2556             : static GEN
    2557         119 : RIpolylog(long m, GEN x, long real, long prec)
    2558             : {
    2559         119 :   GEN y = polylog(m, x, prec);
    2560         119 :   return real? real_i(y): imag_i(y);
    2561             : }
    2562             : GEN
    2563          21 : dilog(GEN x, long prec) { return gpolylog(2, x, prec); }
    2564             : 
    2565             : /* x a floating point number, t_REAL or t_COMPLEX of t_REAL */
    2566             : static GEN
    2567          42 : logabs(GEN x)
    2568             : {
    2569             :   GEN y;
    2570          42 :   if (typ(x) == t_COMPLEX)
    2571             :   {
    2572           7 :     y = logr_abs( cxnorm(x) );
    2573           7 :     shiftr_inplace(y, -1);
    2574             :   } else
    2575          35 :     y = logr_abs(x);
    2576          42 :   return y;
    2577             : }
    2578             : 
    2579             : static GEN
    2580          21 : polylogD(long m, GEN x, long flag, long prec)
    2581             : {
    2582          21 :   long fl = 0, k, l, m2;
    2583             :   pari_sp av;
    2584             :   GEN p1, p2, y;
    2585             : 
    2586          21 :   if (gequal0(x)) return gcopy(x);
    2587          21 :   m2 = m&1;
    2588          21 :   if (gequal1(x) && m>=2) return m2? szeta(m,prec): gen_0;
    2589          21 :   av = avma; l = precision(x);
    2590          21 :   if (!l) { l = prec; x = gtofp(x,l); }
    2591          21 :   p1 = logabs(x);
    2592          21 :   if (signe(p1) > 0) { x = ginv(x); fl = !m2; } else setabssign(p1);
    2593             :   /* |x| <= 1, p1 = - log|x| >= 0 */
    2594          21 :   p2 = gen_1;
    2595          21 :   y = RIpolylog(m, x, m2, l);
    2596          84 :   for (k = 1; k < m; k++)
    2597             :   {
    2598          63 :     GEN t = RIpolylog(m-k, x, m2, l);
    2599          63 :     p2 = gdivgs(gmul(p2,p1), k); /* (-log|x|)^k / k! */
    2600          63 :     y = gadd(y, gmul(p2, t));
    2601             :   }
    2602          21 :   if (m2)
    2603             :   {
    2604          14 :     p1 = flag? gdivgs(p1, -2*m): gdivgs(logabs(gsubsg(1,x)), m);
    2605          14 :     y = gadd(y, gmul(p2, p1));
    2606             :   }
    2607          21 :   if (fl) y = gneg(y);
    2608          21 :   return gerepileupto(av, y);
    2609             : }
    2610             : 
    2611             : static GEN
    2612          14 : polylogP(long m, GEN x, long prec)
    2613             : {
    2614          14 :   long fl = 0, k, l, m2;
    2615             :   pari_sp av;
    2616             :   GEN p1,y;
    2617             : 
    2618          14 :   if (gequal0(x)) return gcopy(x);
    2619          14 :   m2 = m&1;
    2620          14 :   if (gequal1(x) && m>=2) return m2? szeta(m,prec): gen_0;
    2621          14 :   av = avma; l = precision(x);
    2622          14 :   if (!l) { l = prec; x = gtofp(x,l); }
    2623          14 :   p1 = logabs(x);
    2624          14 :   if (signe(p1) > 0) { x = ginv(x); fl = !m2; setsigne(p1, -1); }
    2625             :   /* |x| <= 1 */
    2626          14 :   y = RIpolylog(m, x, m2, l);
    2627          14 :   if (m==1)
    2628             :   {
    2629           7 :     shiftr_inplace(p1, -1); /* log |x| / 2 */
    2630           7 :     y = gadd(y, p1);
    2631             :   }
    2632             :   else
    2633             :   { /* m >= 2, \sum_{0 <= k <= m} 2^k B_k/k! (log |x|)^k Li_{m-k}(x),
    2634             :        with Li_0(x) := -1/2 */
    2635           7 :     GEN u, t = RIpolylog(m-1, x, m2, l);
    2636           7 :     u = gneg_i(p1); /* u = 2 B1 log |x| */
    2637           7 :     y = gadd(y, gmul(u, t));
    2638           7 :     if (m > 2)
    2639             :     {
    2640             :       GEN p2;
    2641           7 :       shiftr_inplace(p1, 1); /* 2log|x| <= 0 */
    2642           7 :       constbern(m>>1);
    2643           7 :       p1 = sqrr(p1);
    2644           7 :       p2 = shiftr(p1,-1);
    2645          21 :       for (k = 2; k < m; k += 2)
    2646             :       {
    2647          14 :         if (k > 2) p2 = divgunu(gmul(p2,p1),k-1); /* 2^k/k! log^k |x|*/
    2648          14 :         t = RIpolylog(m-k, x, m2, l);
    2649          14 :         u = gmul(p2, bernfrac(k));
    2650          14 :         y = gadd(y, gmul(u, t));
    2651             :       }
    2652             :     }
    2653             :   }
    2654          14 :   if (fl) y = gneg(y);
    2655          14 :   return gerepileupto(av, y);
    2656             : }
    2657             : 
    2658             : GEN
    2659         147 : gpolylog(long m, GEN x, long prec)
    2660             : {
    2661         147 :   pari_sp av = avma;
    2662             :   long i, n, v;
    2663             :   GEN a, y;
    2664             : 
    2665         147 :   if (m <= 0)
    2666             :   {
    2667          14 :     a = gmul(x, poleval(eulerianpol(-m, 0), x));
    2668          14 :     return gerepileupto(av, gdiv(a, gpowgs(gsubsg(1, x), 1-m)));
    2669             :   }
    2670         133 :   switch(typ(x))
    2671             :   {
    2672          84 :     case t_INT: case t_REAL: case t_FRAC: case t_COMPLEX: case t_QUAD:
    2673          84 :       return polylog(m,x,prec);
    2674           7 :     case t_POLMOD:
    2675           7 :       return gerepileupto(av, polylogvec(m, polmod_to_embed(x, prec), prec));
    2676           7 :     case t_INTMOD: case t_PADIC: pari_err_IMPL( "padic polylogarithm");
    2677           7 :     case t_VEC: case t_COL: case t_MAT:
    2678           7 :       return polylogvec(m, x, prec);
    2679          28 :     default:
    2680          28 :       av = avma; if (!(y = toser_i(x))) break;
    2681          21 :       if (!m) { set_avma(av); return mkfrac(gen_m1,gen_2); }
    2682          21 :       if (m==1) return gerepileupto(av, Li1(y, prec));
    2683          21 :       if (gequal0(y)) return gerepilecopy(av, y);
    2684          21 :       v = valp(y);
    2685          21 :       if (v < 0) pari_err_DOMAIN("polylog","valuation", "<", gen_0, x);
    2686          14 :       if (v > 0) {
    2687           7 :         n = (lg(y)-3 + v) / v;
    2688           7 :         a = zeroser(varn(y), lg(y)-2);
    2689          35 :         for (i=n; i>=1; i--)
    2690          28 :           a = gmul(y, gadd(a, powis(utoipos(i),-m)));
    2691             :       } else { /* v == 0 */
    2692           7 :         long vy = varn(y);
    2693           7 :         GEN a0 = polcoef(y, 0, -1), t = gdiv(derivser(y), y);
    2694           7 :         a = Li1(y, prec);
    2695          14 :         for (i=2; i<=m; i++)
    2696           7 :           a = gadd(gpolylog(i, a0, prec), integ(gmul(t, a), vy));
    2697             :       }
    2698          14 :       return gerepileupto(av, a);
    2699             :   }
    2700           7 :   pari_err_TYPE("gpolylog",x);
    2701             :   return NULL; /* LCOV_EXCL_LINE */
    2702             : }
    2703             : 
    2704             : GEN
    2705         133 : polylog0(long m, GEN x, long flag, long prec)
    2706             : {
    2707         133 :   switch(flag)
    2708             :   {
    2709          91 :     case 0: return gpolylog(m,x,prec);
    2710          14 :     case 1: return polylogD(m,x,0,prec);
    2711           7 :     case 2: return polylogD(m,x,1,prec);
    2712          14 :     case 3: return polylogP(m,x,prec);
    2713           7 :     default: pari_err_FLAG("polylog");
    2714             :   }
    2715             :   return NULL; /* LCOV_EXCL_LINE */
    2716             : }
    2717             : 
    2718             : GEN
    2719       22379 : upper_to_cx(GEN x, long *prec)
    2720             : {
    2721       22379 :   long tx = typ(x), l;
    2722       22379 :   if (tx == t_QUAD) { x = quadtofp(x, *prec); tx = typ(x); }
    2723       22379 :   switch(tx)
    2724             :   {
    2725       22358 :     case t_COMPLEX:
    2726       22358 :       if (gsigne(gel(x,2)) > 0) break; /*fall through*/
    2727             :     case t_REAL: case t_INT: case t_FRAC:
    2728          14 :       pari_err_DOMAIN("modular function", "Im(argument)", "<=", gen_0, x);
    2729           7 :     default:
    2730           7 :       pari_err_TYPE("modular function", x);
    2731             :   }
    2732       22358 :   l = precision(x); if (l) *prec = l;
    2733       22358 :   return x;
    2734             : }
    2735             : 
    2736             : /* sqrt(3)/2 */
    2737             : static GEN
    2738        1456 : sqrt32(long prec) { GEN z = sqrtr_abs(utor(3,prec)); setexpo(z, -1); return z; }
    2739             : /* exp(i x), x = k pi/12 */
    2740             : static GEN
    2741        2392 : e12(ulong k, long prec)
    2742             : {
    2743             :   int s, sPi, sPiov2;
    2744             :   GEN z, t;
    2745        2392 :   k %= 24;
    2746        2392 :   if (!k) return gen_1;
    2747        2392 :   if (k == 12) return gen_m1;
    2748        2392 :   if (k >12) { s = 1; k = 24 - k; } else s = 0; /* x -> 2pi - x */
    2749        2392 :   if (k > 6) { sPi = 1; k = 12 - k; } else sPi = 0; /* x -> pi  - x */
    2750        2392 :   if (k > 3) { sPiov2 = 1; k = 6 - k; } else sPiov2 = 0; /* x -> pi/2 - x */
    2751        2392 :   z = cgetg(3, t_COMPLEX);
    2752        2392 :   switch(k)
    2753             :   {
    2754         439 :     case 0: gel(z,1) = icopy(gen_1); gel(z,2) = gen_0; break;
    2755         574 :     case 1: t = gmul2n(addrs(sqrt32(prec), 1), -1);
    2756         574 :       gel(z,1) = sqrtr(t);
    2757         574 :       gel(z,2) = gmul2n(invr(gel(z,1)), -2); break;
    2758             : 
    2759         882 :     case 2: gel(z,1) = sqrt32(prec);
    2760         882 :             gel(z,2) = real2n(-1, prec); break;
    2761             : 
    2762         497 :     case 3: gel(z,1) = sqrtr_abs(real2n(-1,prec));
    2763         497 :             gel(z,2) = rcopy(gel(z,1)); break;
    2764             :   }
    2765        2392 :   if (sPiov2) swap(gel(z,1), gel(z,2));
    2766        2392 :   if (sPi) togglesign(gel(z,1));
    2767        2392 :   if (s)   togglesign(gel(z,2));
    2768        2392 :   return z;
    2769             : }
    2770             : /* z a t_FRAC */
    2771             : static GEN
    2772       11058 : expIPifrac(GEN z, long prec)
    2773             : {
    2774       11058 :   GEN n = gel(z,1), d = gel(z,2);
    2775       11058 :   ulong r, q = uabsdivui_rem(12, d, &r);
    2776       11058 :   if (!r) return e12(q * umodiu(n, 24), prec); /* 12 | d */
    2777        8666 :   n = centermodii(n, shifti(d,1), d);
    2778        8666 :   return expIr(divri(mulri(mppi(prec), n), d));
    2779             : }
    2780             : /* exp(i Pi z), z a t_INT or t_FRAC */
    2781             : static GEN
    2782       10997 : expIPiQ(GEN z, long prec)
    2783             : {
    2784       10997 :   if (typ(z) == t_INT) return mpodd(z)? gen_m1: gen_1;
    2785        1946 :   return expIPifrac(z, prec);
    2786             : }
    2787             : 
    2788             : /* convert power of 2 t_REAL to rational */
    2789             : static GEN
    2790        1132 : real2nQ(GEN x)
    2791             : {
    2792        1132 :   long e = expo(x);
    2793             :   GEN z;
    2794        1132 :   if (e < 0)
    2795         474 :     z = mkfrac(signe(x) < 0? gen_m1: gen_1, int2n(-e));
    2796             :   else
    2797             :   {
    2798         658 :     z = int2n(e);
    2799         658 :     if (signe(x) < 0) togglesign_safe(&z);
    2800             :   }
    2801        1132 :   return z;
    2802             : }
    2803             : /* x a real number */
    2804             : GEN
    2805         656 : expIPiR(GEN x, long prec)
    2806             : {
    2807         656 :   if (typ(x) == t_REAL && absrnz_equal2n(x)) x = real2nQ(x);
    2808         656 :   switch(typ(x))
    2809             :   {
    2810           0 :     case t_INT:  return mpodd(x)? gen_m1: gen_1;
    2811         452 :     case t_FRAC: return expIPifrac(x, prec);
    2812             :   }
    2813         204 :   return expIr(mulrr(mppi(prec), x));
    2814             : }
    2815             : /* z a t_COMPLEX */
    2816             : GEN
    2817       50701 : expIPiC(GEN z, long prec)
    2818             : {
    2819             :   GEN pi, r, x, y;
    2820       50701 :   if (typ(z) != t_COMPLEX) return expIPiR(z, prec);
    2821       50351 :   x = gel(z,1);
    2822       50351 :   y = gel(z,2); if (gequal0(y)) return expIPiR(x, prec);
    2823       50045 :   pi = mppi(prec);
    2824       50045 :   r = gmul(pi, y); togglesign(r); r = mpexp(r); /* exp(-pi y) */
    2825       50045 :   if (typ(x) == t_REAL && absrnz_equal2n(x)) x = real2nQ(x);
    2826       50045 :   switch(typ(x))
    2827             :   {
    2828        2702 :     case t_INT: if (mpodd(x)) togglesign(r);
    2829        2702 :                 return r;
    2830        8660 :     case t_FRAC: return gmul(r, expIPifrac(x, prec));
    2831             :   }
    2832       38683 :   return gmul(r, expIr(mulrr(pi, x)));
    2833             : }
    2834             : /* exp(I x y), more efficient for x in R, y pure imaginary */
    2835             : GEN
    2836         133 : expIxy(GEN x, GEN y, long prec) { return gexp(gmul(x, mulcxI(y)), prec); }
    2837             : 
    2838             : static GEN
    2839       13517 : qq(GEN x, long prec)
    2840             : {
    2841       13517 :   long tx = typ(x);
    2842             :   GEN y;
    2843             : 
    2844       13517 :   if (is_scalar_t(tx))
    2845             :   {
    2846       13475 :     if (tx == t_PADIC) return x;
    2847       13461 :     x = upper_to_cx(x, &prec);
    2848       13447 :     return cxtoreal(expIPiC(gmul2n(x,1), prec)); /* e(x) */
    2849             :   }
    2850          42 :   if (! ( y = toser_i(x)) ) pari_err_TYPE("modular function", x);
    2851          42 :   return y;
    2852             : }
    2853             : 
    2854             : /* return (y * X^d) + x. Assume d > 0, x != 0, valp(x) = 0 */
    2855             : static GEN
    2856          21 : ser_addmulXn(GEN y, GEN x, long d)
    2857             : {
    2858          21 :   long i, lx, ly, l = valp(y) + d; /* > 0 */
    2859             :   GEN z;
    2860             : 
    2861          21 :   lx = lg(x);
    2862          21 :   ly = lg(y) + l; if (lx < ly) ly = lx;
    2863          21 :   if (l > lx-2) return gcopy(x);
    2864          21 :   z = cgetg(ly,t_SER);
    2865          77 :   for (i=2; i<=l+1; i++) gel(z,i) = gel(x,i);
    2866          70 :   for (   ; i < ly; i++) gel(z,i) = gadd(gel(x,i),gel(y,i-l));
    2867          21 :   z[1] = x[1]; return z;
    2868             : }
    2869             : 
    2870             : /* q a t_POL s.t. q(0) != 0, v > 0, Q = x^v*q; return \prod_i (1-Q^i) */
    2871             : static GEN
    2872          28 : RgXn_eta(GEN q, long v, long lim)
    2873             : {
    2874          28 :   pari_sp av = avma;
    2875             :   GEN qn, ps, y;
    2876             :   ulong vps, vqn, n;
    2877             : 
    2878          28 :   if (!degpol(q) && isint1(gel(q,2))) return eta_ZXn(v, lim+v);
    2879           7 :   y = qn = ps = pol_1(0);
    2880           7 :   vps = vqn = 0;
    2881           7 :   for(n = 0;; n++)
    2882           7 :   { /* qn = q^n,  ps = (-1)^n q^(n(3n+1)/2),
    2883             :      * vps, vqn valuation of ps, qn HERE */
    2884          14 :     pari_sp av2 = avma;
    2885          14 :     ulong vt = vps + 2*vqn + v; /* valuation of t at END of loop body */
    2886             :     long k1, k2;
    2887             :     GEN t;
    2888          14 :     vqn += v; vps = vt + vqn; /* valuation of qn, ps at END of body */
    2889          14 :     k1 = lim + v - vt + 1;
    2890          14 :     k2 = k1 - vqn; /* = lim + v - vps + 1 */
    2891          14 :     if (k1 <= 0) break;
    2892          14 :     t = RgX_mul(q, RgX_sqr(qn));
    2893          14 :     t = RgXn_red_shallow(t, k1);
    2894          14 :     t = RgX_mul(ps,t);
    2895          14 :     t = RgXn_red_shallow(t, k1);
    2896          14 :     t = RgX_neg(t); /* t = (-1)^(n+1) q^(n(3n+1)/2 + 2n+1) */
    2897          14 :     t = gerepileupto(av2, t);
    2898          14 :     y = RgX_addmulXn_shallow(t, y, vt);
    2899          14 :     if (k2 <= 0) break;
    2900             : 
    2901           7 :     qn = RgX_mul(qn,q);
    2902           7 :     ps = RgX_mul(t,qn);
    2903           7 :     ps = RgXn_red_shallow(ps, k2);
    2904           7 :     y = RgX_addmulXn_shallow(ps, y, vps);
    2905             : 
    2906           7 :     if (gc_needed(av,1))
    2907             :     {
    2908           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"eta, n = %ld", n);
    2909           0 :       gerepileall(av, 3, &y, &qn, &ps);
    2910             :     }
    2911             :   }
    2912           7 :   return y;
    2913             : }
    2914             : 
    2915             : static GEN
    2916       16443 : inteta(GEN q)
    2917             : {
    2918       16443 :   long tx = typ(q);
    2919             :   GEN ps, qn, y;
    2920             : 
    2921       16443 :   y = gen_1; qn = gen_1; ps = gen_1;
    2922       16443 :   if (tx==t_PADIC)
    2923             :   {
    2924          28 :     if (valp(q) <= 0) pari_err_DOMAIN("eta", "v_p(q)", "<=",gen_0,q);
    2925             :     for(;;)
    2926          56 :     {
    2927          77 :       GEN t = gneg_i(gmul(ps,gmul(q,gsqr(qn))));
    2928          77 :       y = gadd(y,t); qn = gmul(qn,q); ps = gmul(t,qn);
    2929          77 :       t = y;
    2930          77 :       y = gadd(y,ps); if (gequal(t,y)) return y;
    2931             :     }
    2932             :   }
    2933             : 
    2934       16415 :   if (tx == t_SER)
    2935             :   {
    2936             :     ulong vps, vqn;
    2937          42 :     long l = lg(q), v, n;
    2938             :     pari_sp av;
    2939             : 
    2940          42 :     v = valp(q); /* handle valuation separately to avoid overflow */
    2941          42 :     if (v <= 0) pari_err_DOMAIN("eta", "v_p(q)", "<=",gen_0,q);
    2942          35 :     y = ser2pol_i(q, l); /* t_SER inefficient when input has low degree */
    2943          35 :     n = degpol(y);
    2944          35 :     if (n <= (l>>2))
    2945             :     {
    2946          28 :       GEN z = RgXn_eta(y, v, l-2);
    2947          28 :       setvarn(z, varn(y)); return RgX_to_ser(z, l+v);
    2948             :     }
    2949           7 :     q = leafcopy(q); av = avma;
    2950           7 :     setvalp(q, 0);
    2951           7 :     y = scalarser(gen_1, varn(q), l+v);
    2952           7 :     vps = vqn = 0;
    2953           7 :     for(n = 0;; n++)
    2954           7 :     { /* qn = q^n,  ps = (-1)^n q^(n(3n+1)/2) */
    2955          14 :       ulong vt = vps + 2*vqn + v;
    2956             :       long k;
    2957             :       GEN t;
    2958          14 :       t = gneg_i(gmul(ps,gmul(q,gsqr(qn))));
    2959             :       /* t = (-1)^(n+1) q^(n(3n+1)/2 + 2n+1) */
    2960          14 :       y = ser_addmulXn(t, y, vt);
    2961          14 :       vqn += v; vps = vt + vqn;
    2962          14 :       k = l+v - vps; if (k <= 2) return y;
    2963             : 
    2964           7 :       qn = gmul(qn,q); ps = gmul(t,qn);
    2965           7 :       y = ser_addmulXn(ps, y, vps);
    2966           7 :       setlg(q, k);
    2967           7 :       setlg(qn, k);
    2968           7 :       setlg(ps, k);
    2969           7 :       if (gc_needed(av,3))
    2970             :       {
    2971           0 :         if(DEBUGMEM>1) pari_warn(warnmem,"eta");
    2972           0 :         gerepileall(av, 3, &y, &qn, &ps);
    2973             :       }
    2974             :     }
    2975             :   }
    2976             :   {
    2977       16373 :     long l = -prec2nbits(precision(q));
    2978       16373 :     pari_sp av = avma;
    2979             : 
    2980             :     for(;;)
    2981       52899 :     {
    2982       69272 :       GEN t = gneg_i(gmul(ps,gmul(q,gsqr(qn))));
    2983             :       /* qn = q^n
    2984             :        * ps = (-1)^n q^(n(3n+1)/2)
    2985             :        * t = (-1)^(n+1) q^(n(3n+1)/2 + 2n+1) */
    2986       69272 :       y = gadd(y,t); qn = gmul(qn,q); ps = gmul(t,qn);
    2987       69272 :       y = gadd(y,ps);
    2988       69272 :       if (gexpo(ps)-gexpo(y) < l) return y;
    2989       52899 :       if (gc_needed(av,3))
    2990             :       {
    2991           0 :         if(DEBUGMEM>1) pari_warn(warnmem,"eta");
    2992           0 :         gerepileall(av, 3, &y, &qn, &ps);
    2993             :       }
    2994             :     }
    2995             :   }
    2996             : }
    2997             : 
    2998             : GEN
    2999          77 : eta(GEN x, long prec)
    3000             : {
    3001          77 :   pari_sp av = avma;
    3002          77 :   GEN z = inteta( qq(x,prec) );
    3003          49 :   if (typ(z) == t_SER) return gerepilecopy(av, z);
    3004          14 :   return gerepileupto(av, z);
    3005             : }
    3006             : 
    3007             : /* s(h,k) = sum(n = 1, k-1, (n/k)*(frac(h*n/k) - 1/2))
    3008             :  * Knuth's algorithm. h integer, k integer > 0, (h,k) = 1 */
    3009             : GEN
    3010        6258 : sumdedekind_coprime(GEN h, GEN k)
    3011             : {
    3012        6258 :   pari_sp av = avma;
    3013             :   GEN s2, s1, p, pp;
    3014             :   long s;
    3015        6258 :   if (lgefint(k) == 3 && uel(k,2) <= (2*(ulong)LONG_MAX) / 3)
    3016             :   {
    3017        6251 :     ulong kk = k[2], hh = umodiu(h, kk);
    3018             :     long s1, s2;
    3019             :     GEN v;
    3020        6251 :     if (signe(k) < 0) { k = negi(k); hh = Fl_neg(hh, kk); }
    3021        6251 :     v = u_sumdedekind_coprime(hh, kk);
    3022        6251 :     s1 = v[1]; s2 = v[2];
    3023        6251 :     return gerepileupto(av, gdiv(addis(mulis(k,s1), s2), muluu(12, kk)));
    3024             :   }
    3025           7 :   s = 1;
    3026           7 :   s1 = gen_0; p = gen_1; pp = gen_0;
    3027           7 :   s2 = h = modii(h, k);
    3028          35 :   while (signe(h)) {
    3029          28 :     GEN r, nexth, a = dvmdii(k, h, &nexth);
    3030          28 :     if (is_pm1(h)) s2 = s == 1? addii(s2, p): subii(s2, p);
    3031          28 :     s1 = s == 1? addii(s1, a): subii(s1, a);
    3032          28 :     s = -s;
    3033          28 :     k = h; h = nexth;
    3034          28 :     r = addii(mulii(a,p), pp); pp = p; p = r;
    3035             :   }
    3036             :   /* at this point p = original k */
    3037           7 :   if (s == -1) s1 = subiu(s1, 3);
    3038           7 :   return gerepileupto(av, gdiv(addii(mulii(p,s1), s2), muliu(p,12)));
    3039             : }
    3040             : /* as above, for ulong arguments.
    3041             :  * k integer > 0, 0 <= h < k, (h,k) = 1. Returns [s1,s2] such that
    3042             :  * s(h,k) = (s2 + k s1) / (12k). Requires max(h + k/2, k) < LONG_MAX
    3043             :  * to avoid overflow, in particular k <= LONG_MAX * 2/3 is fine */
    3044             : GEN
    3045        6251 : u_sumdedekind_coprime(long h, long k)
    3046             : {
    3047        6251 :   long s = 1, s1 = 0, s2 = h, p = 1, pp = 0;
    3048       11424 :   while (h) {
    3049        5173 :     long r, nexth = k % h, a = k / h; /* a >= 1, a >= 2 if h = 1 */
    3050        5173 :     if (h == 1) s2 += p * s; /* occurs exactly once, last step */
    3051        5173 :     s1 += a * s;
    3052        5173 :     s = -s;
    3053        5173 :     k = h; h = nexth;
    3054        5173 :     r = a*p + pp; pp = p; p = r; /* p >= pp >= 0 */
    3055             :   }
    3056             :   /* in the above computation, p increases from 1 to original k,
    3057             :    * -k/2 <= s2 <= h + k/2, and |s1| <= k */
    3058        6251 :   if (s < 0) s1 -= 3; /* |s1| <= k+3 ? */
    3059             :   /* But in fact, |s2 + p s1| <= k^2 + 1/2 - 3k; if (s < 0), we have
    3060             :    * |s2| <= k/2 and it follows that |s1| < k here as well */
    3061             :   /* p = k; s(h,k) = (s2 + p s1)/12p. */
    3062        6251 :   return mkvecsmall2(s1, s2);
    3063             : }
    3064             : GEN
    3065          28 : sumdedekind(GEN h, GEN k)
    3066             : {
    3067          28 :   pari_sp av = avma;
    3068             :   GEN d;
    3069          28 :   if (typ(h) != t_INT) pari_err_TYPE("sumdedekind",h);
    3070          28 :   if (typ(k) != t_INT) pari_err_TYPE("sumdedekind",k);
    3071          28 :   d = gcdii(h,k);
    3072          28 :   if (!is_pm1(d))
    3073             :   {
    3074           7 :     h = diviiexact(h, d);
    3075           7 :     k = diviiexact(k, d);
    3076             :   }
    3077          28 :   return gerepileupto(av, sumdedekind_coprime(h,k));
    3078             : }
    3079             : 
    3080             : /* eta(x); assume Im x >> 0 (e.g. x in SL2's standard fundamental domain) */
    3081             : static GEN
    3082       17297 : eta_reduced(GEN x, long prec)
    3083             : {
    3084       17297 :   GEN z = expIPiC(gdivgs(x, 12), prec); /* e(x/24) */
    3085       17297 :   if (24 * gexpo(z) >= -prec2nbits(prec))
    3086       16352 :     z = gmul(z, inteta( gpowgs(z,24) ));
    3087       17297 :   return z;
    3088             : }
    3089             : 
    3090             : /* x = U.z (flag = 1), or x = U^(-1).z (flag = 0)
    3091             :  * Return [s,t] such that eta(z) = eta(x) * sqrt(s) * exp(I Pi t) */
    3092             : static GEN
    3093       17311 : eta_correction(GEN x, GEN U, long flag)
    3094             : {
    3095             :   GEN a,b,c,d, s,t;
    3096             :   long sc;
    3097       17311 :   a = gcoeff(U,1,1);
    3098       17311 :   b = gcoeff(U,1,2);
    3099       17311 :   c = gcoeff(U,2,1);
    3100       17311 :   d = gcoeff(U,2,2);
    3101             :   /* replace U by U^(-1) */
    3102       17311 :   if (flag) {
    3103        8939 :     swap(a,d);
    3104        8939 :     togglesign_safe(&b);
    3105        8939 :     togglesign_safe(&c);
    3106             :   }
    3107       17311 :   sc = signe(c);
    3108       17311 :   if (!sc) {
    3109       11081 :     if (signe(d) < 0) togglesign_safe(&b);
    3110       11081 :     s = gen_1;
    3111       11081 :     t = sstoQ(umodiu(b, 24), 12);
    3112             :   } else {
    3113        6230 :     if (sc < 0) {
    3114        1722 :       togglesign_safe(&a);
    3115        1722 :       togglesign_safe(&b);
    3116        1722 :       togglesign_safe(&c);
    3117        1722 :       togglesign_safe(&d);
    3118             :     } /* now c > 0 */
    3119        6230 :     s = mulcxmI(gadd(gmul(c,x), d));
    3120        6230 :     t = gadd(gdiv(addii(a,d),muliu(c,12)), sumdedekind_coprime(negi(d),c));
    3121             :     /* correction : exp(I Pi (((a+d)/12c) + s(-d,c)) ) sqrt(-i(cx+d))  */
    3122             :   }
    3123       17311 :   return mkvec2(s, t);
    3124             : }
    3125             : 
    3126             : /* returns the true value of eta(x) for Im(x) > 0, using reduction to
    3127             :  * standard fundamental domain */
    3128             : GEN
    3129        8869 : trueeta(GEN x, long prec)
    3130             : {
    3131        8869 :   pari_sp av = avma;
    3132             :   GEN U, st, s, t;
    3133             : 
    3134        8869 :   if (!is_scalar_t(typ(x))) pari_err_TYPE("trueeta",x);
    3135        8869 :   x = upper_to_cx(x, &prec);
    3136        8869 :   x = cxredsl2(x, &U);
    3137        8869 :   st = eta_correction(x, U, 1);
    3138        8869 :   x = eta_reduced(x, prec);
    3139        8869 :   s = gel(st, 1);
    3140        8869 :   t = gel(st, 2);
    3141        8869 :   x = gmul(x, expIPiQ(t, prec));
    3142        8869 :   if (s != gen_1) x = gmul(x, gsqrt(s, prec));
    3143        8869 :   return gerepileupto(av, x);
    3144             : }
    3145             : 
    3146             : GEN
    3147          77 : eta0(GEN x, long flag,long prec)
    3148          77 : { return flag? trueeta(x,prec): eta(x,prec); }
    3149             : 
    3150             : /* eta(q) = 1 + \sum_{n>0} (-1)^n * (q^(n(3n-1)/2) + q^(n(3n+1)/2)) */
    3151             : static GEN
    3152           7 : ser_eta(long prec)
    3153             : {
    3154           7 :   GEN e = cgetg(prec+2, t_SER), ed = e+2;
    3155             :   long n, j;
    3156           7 :   e[1] = evalsigne(1)|_evalvalp(0)|evalvarn(0);
    3157           7 :   gel(ed,0) = gen_1;
    3158         483 :   for (n = 1; n < prec; n++) gel(ed,n) = gen_0;
    3159          49 :   for (n = 1, j = 0; n < prec; n++)
    3160             :   {
    3161             :     GEN s;
    3162          49 :     j += 3*n-2; /* = n*(3*n-1) / 2 */;
    3163          49 :     if (j >= prec) break;
    3164          42 :     s = odd(n)? gen_m1: gen_1;
    3165          42 :     gel(ed, j) = s;
    3166          42 :     if (j+n >= prec) break;
    3167          42 :     gel(ed, j+n) = s;
    3168             :   }
    3169           7 :   return e;
    3170             : }
    3171             : 
    3172             : static GEN
    3173         476 : coeffEu(GEN fa)
    3174             : {
    3175         476 :   pari_sp av = avma;
    3176         476 :   return gerepileuptoint(av, mului(65520, usumdivk_fact(fa,11)));
    3177             : }
    3178             : /* E12 = 1 + q*E/691 */
    3179             : static GEN
    3180           7 : ser_E(long prec)
    3181             : {
    3182           7 :   GEN e = cgetg(prec+2, t_SER), ed = e+2;
    3183           7 :   GEN F = vecfactoru_i(2, prec); /* F[n] = factoru(n+1) */
    3184             :   long n;
    3185           7 :   e[1] = evalsigne(1)|_evalvalp(0)|evalvarn(0);
    3186           7 :   gel(ed,0) = utoipos(65520);
    3187         483 :   for (n = 1; n < prec; n++) gel(ed,n) = coeffEu(gel(F,n));
    3188           7 :   return e;
    3189             : }
    3190             : /* j = E12/Delta + 432000/691, E12 = 1 + q*E/691 */
    3191             : static GEN
    3192           7 : ser_j2(long prec, long v)
    3193             : {
    3194           7 :   pari_sp av = avma;
    3195           7 :   GEN iD = gpowgs(ginv(ser_eta(prec)), 24); /* q/Delta */
    3196           7 :   GEN J = gmul(ser_E(prec), iD);
    3197           7 :   setvalp(iD,-1); /* now 1/Delta */
    3198           7 :   J = gadd(gdivgs(J, 691), iD);
    3199           7 :   J = gerepileupto(av, J);
    3200           7 :   if (prec > 1) gel(J,3) = utoipos(744);
    3201           7 :   setvarn(J,v); return J;
    3202             : }
    3203             : 
    3204             : /* j(q) = \sum_{n >= -1} c(n)q^n,
    3205             :  * \sum_{n = -1}^{N-1} c(n) (-10n \sigma_3(N-n) + 21 \sigma_5(N-n))
    3206             :  * = c(N) (N+1)/24 */
    3207             : static GEN
    3208          14 : ser_j(long prec, long v)
    3209             : {
    3210             :   GEN j, J, S3, S5, F;
    3211             :   long i, n;
    3212          14 :   if (prec > 64) return ser_j2(prec, v);
    3213           7 :   S3 = cgetg(prec+1, t_VEC);
    3214           7 :   S5 = cgetg(prec+1,t_VEC);
    3215           7 :   F = vecfactoru_i(1, prec);
    3216          35 :   for (n = 1; n <= prec; n++)
    3217             :   {
    3218          28 :     GEN fa = gel(F,n);
    3219          28 :     gel(S3,n) = mului(10, usumdivk_fact(fa,3));
    3220          28 :     gel(S5,n) = mului(21, usumdivk_fact(fa,5));
    3221             :   }
    3222           7 :   J = cgetg(prec+2, t_SER),
    3223           7 :   J[1] = evalvarn(v)|evalsigne(1)|evalvalp(-1);
    3224           7 :   j = J+3;
    3225           7 :   gel(j,-1) = gen_1;
    3226           7 :   gel(j,0) = utoipos(744);
    3227           7 :   gel(j,1) = utoipos(196884);
    3228          21 :   for (n = 2; n < prec; n++)
    3229             :   {
    3230          14 :     pari_sp av = avma;
    3231          14 :     GEN c, s3 = gel(S3,n+1), s5 = gel(S5,n+1);
    3232          14 :     c = addii(s3, s5);
    3233          49 :     for (i = 0; i < n; i++)
    3234             :     {
    3235          35 :       s3 = gel(S3,n-i); s5 = gel(S5,n-i);
    3236          35 :       c = addii(c, mulii(gel(j,i), subii(s5, mului(i,s3))));
    3237             :     }
    3238          14 :     gel(j,n) = gerepileuptoint(av, diviuexact(muliu(c,24), n+1));
    3239             :   }
    3240           7 :   return J;
    3241             : }
    3242             : 
    3243             : GEN
    3244          42 : jell(GEN x, long prec)
    3245             : {
    3246          42 :   long tx = typ(x);
    3247          42 :   pari_sp av = avma;
    3248             :   GEN q, h, U;
    3249             : 
    3250          42 :   if (!is_scalar_t(tx))
    3251             :   {
    3252             :     long v;
    3253          21 :     if (gequalX(x)) return ser_j(precdl, varn(x));
    3254          21 :     q = toser_i(x); if (!q) pari_err_TYPE("ellj",x);
    3255          14 :     v = fetch_var_higher();
    3256          14 :     h = ser_j(lg(q)-2, v);
    3257          14 :     h = gsubst(h, v, q);
    3258          14 :     delete_var(); return gerepileupto(av, h);
    3259             :   }
    3260          21 :   if (tx == t_PADIC)
    3261             :   {
    3262           7 :     GEN p2, p1 = gdiv(inteta(gsqr(x)), inteta(x));
    3263           7 :     p1 = gmul2n(gsqr(p1),1);
    3264           7 :     p1 = gmul(x,gpowgs(p1,12));
    3265           7 :     p2 = gaddsg(768,gadd(gsqr(p1),gdivsg(4096,p1)));
    3266           7 :     p1 = gmulsg(48,p1);
    3267           7 :     return gerepileupto(av, gadd(p2,p1));
    3268             :   }
    3269             :   /* Let h = Delta(2x) / Delta(x), then j(x) = (1 + 256h)^3 / h */
    3270          14 :   x = upper_to_cx(x, &prec);
    3271           7 :   x = cxredsl2(x, &U); /* forget about Ua : j has weight 0 */
    3272             :   { /* cf eta_reduced, raised to power 24
    3273             :      * Compute
    3274             :      *   t = (inteta(q(2x)) / inteta(q(x))) ^ 24;
    3275             :      * then
    3276             :      *   h = t * (q(2x) / q(x) = t * q(x);
    3277             :      * but inteta(q) costly and useless if expo(q) << 1  => inteta(q) = 1.
    3278             :      * log_2 ( exp(-2Pi Im tau) ) < -prec2nbits(prec)
    3279             :      * <=> Im tau > prec2nbits(prec) * log(2) / 2Pi */
    3280           7 :     long C = (long)prec2nbits_mul(prec, M_LN2/(2*M_PI));
    3281           7 :     q = expIPiC(gmul2n(x,1), prec); /* e(x) */
    3282           7 :     if (gcmpgs(gel(x,2), C) > 0) /* eta(q(x)) = 1 : no need to compute q(2x) */
    3283           0 :       h = q;
    3284             :     else
    3285             :     {
    3286           7 :       GEN t = gdiv(inteta(gsqr(q)), inteta(q));
    3287           7 :       h = gmul(q, gpowgs(t, 24));
    3288             :     }
    3289             :   }
    3290             :   /* real_1 important ! gaddgs(, 1) could increase the accuracy ! */
    3291           7 :   return gerepileupto(av, gdiv(gpowgs(gadd(gmul2n(h,8), real_1(prec)), 3), h));
    3292             : }
    3293             : 
    3294             : static GEN
    3295        8372 : to_form(GEN a, GEN w, GEN C)
    3296        8372 : { return mkvec3(a, w, diviiexact(C, a)); }
    3297             : static GEN
    3298        8372 : form_to_quad(GEN f, GEN sqrtD)
    3299             : {
    3300        8372 :   long a = itos(gel(f,1)), a2 = a << 1;
    3301        8372 :   GEN b = gel(f,2);
    3302        8372 :   return mkcomplex(gdivgs(b, -a2), gdivgs(sqrtD, a2));
    3303             : }
    3304             : static GEN
    3305        8372 : eta_form(GEN f, GEN sqrtD, GEN *s_t, long prec)
    3306             : {
    3307        8372 :   GEN U, t = form_to_quad(redimagsl2(f, &U), sqrtD);
    3308        8372 :   *s_t = eta_correction(t, U, 0);
    3309        8372 :   return eta_reduced(t, prec);
    3310             : }
    3311             : 
    3312             : /* eta(t/p)eta(t/q) / (eta(t)eta(t/pq)), t = (-w + sqrt(D)) / 2a */
    3313             : GEN
    3314        2093 : double_eta_quotient(GEN a, GEN w, GEN D, long p, long q, GEN pq, GEN sqrtD)
    3315             : {
    3316        2093 :   GEN C = shifti(subii(sqri(w), D), -2);
    3317             :   GEN d, t, z, zp, zq, zpq, s_t, s_tp, s_tpq, s, sp, spq;
    3318        2093 :   long prec = realprec(sqrtD);
    3319             : 
    3320        2093 :   z = eta_form(to_form(a, w, C), sqrtD, &s_t, prec);
    3321        2093 :   s = gel(s_t, 1);
    3322        2093 :   zp = eta_form(to_form(mului(p, a), w, C), sqrtD, &s_tp, prec);
    3323        2093 :   sp = gel(s_tp, 1);
    3324        2093 :   zpq = eta_form(to_form(mulii(pq, a), w, C), sqrtD, &s_tpq, prec);
    3325        2093 :   spq = gel(s_tpq, 1);
    3326        2093 :   if (p == q) {
    3327           0 :     z = gdiv(gsqr(zp), gmul(z, zpq));
    3328           0 :     t = gsub(gmul2n(gel(s_tp,2), 1),
    3329           0 :              gadd(gel(s_t,2), gel(s_tpq,2)));
    3330           0 :     if (sp != gen_1) z = gmul(z, sp);
    3331             :   } else {
    3332             :     GEN s_tq, sq;
    3333        2093 :     zq = eta_form(to_form(mului(q, a), w, C), sqrtD, &s_tq, prec);
    3334        2093 :     sq = gel(s_tq, 1);
    3335        2093 :     z = gdiv(gmul(zp, zq), gmul(z, zpq));
    3336        2093 :     t = gsub(gadd(gel(s_tp,2), gel(s_tq,2)),
    3337        2093 :              gadd(gel(s_t,2), gel(s_tpq,2)));
    3338        2093 :     if (sp != gen_1) z = gmul(z, gsqrt(sp, prec));
    3339        2093 :     if (sq != gen_1) z = gmul(z, gsqrt(sq, prec));
    3340             :   }
    3341        2093 :   d = NULL;
    3342        2093 :   if (s != gen_1) d = gsqrt(s, prec);
    3343        2093 :   if (spq != gen_1) {
    3344        2065 :     GEN x = gsqrt(spq, prec);
    3345        2065 :     d = d? gmul(d, x): x;
    3346             :   }
    3347        2093 :   if (d) z = gdiv(z, d);
    3348        2093 :   return gmul(z, expIPiQ(t, prec));
    3349             : }
    3350             : 
    3351             : typedef struct { GEN u; long v, t; } cxanalyze_t;
    3352             : 
    3353             : /* Check whether a t_COMPLEX, t_REAL or t_INT z != 0 can be written as
    3354             :  * z = u * 2^(v/2) * exp(I Pi/4 t), u > 0, v = 0,1 and -3 <= t <= 4.
    3355             :  * Allow z t_INT/t_REAL to simplify handling of eta_correction() output */
    3356             : static int
    3357          70 : cxanalyze(cxanalyze_t *T, GEN z)
    3358             : {
    3359             :   GEN a, b;
    3360             :   long ta, tb;
    3361             : 
    3362          70 :   T->v = 0;
    3363          70 :   if (is_intreal_t(typ(z)))
    3364             :   {
    3365          63 :     T->u = mpabs_shallow(z);
    3366          63 :     T->t = signe(z) < 0? 4: 0;
    3367          63 :     return 1;
    3368             :   }
    3369           7 :   a = gel(z,1); ta = typ(a);
    3370           7 :   b = gel(z,2); tb = typ(b);
    3371             : 
    3372           7 :   T->t = 0;
    3373           7 :   if (ta == t_INT && !signe(a))
    3374             :   {
    3375           0 :     T->u = R_abs_shallow(b);
    3376           0 :     T->t = gsigne(b) < 0? -2: 2;
    3377           0 :     return 1;
    3378             :   }
    3379           7 :   if (tb == t_INT && !signe(b))
    3380             :   {
    3381           0 :     T->u = R_abs_shallow(a);
    3382           0 :     T->t = gsigne(a) < 0? 4: 0;
    3383           0 :     return 1;
    3384             :   }
    3385           7 :   if (ta != tb || ta == t_REAL) { T->u = z; return 0; }
    3386             :   /* a,b both non zero, both t_INT or t_FRAC */
    3387           7 :   if (ta == t_INT)
    3388             :   {
    3389           7 :     if (!absequalii(a, b)) return 0;
    3390           7 :     T->u = absi_shallow(a);
    3391           7 :     T->v = 1;
    3392           7 :     if (signe(a) == signe(b))
    3393           0 :     { T->t = signe(a) < 0? -3: 1; }
    3394             :     else
    3395           7 :     { T->t = signe(a) < 0? 3: -1; }
    3396             :   }
    3397             :   else
    3398             :   {
    3399           0 :     if (!absequalii(gel(a,2), gel(b,2)) || !absequalii(gel(a,1),gel(b,1)))
    3400           0 :       return 0;
    3401           0 :     T->u = absfrac_shallow(a);
    3402           0 :     T->v = 1;
    3403           0 :     a = gel(a,1);
    3404           0 :     b = gel(b,1);
    3405           0 :     if (signe(a) == signe(b))
    3406           0 :     { T->t = signe(a) < 0? -3: 1; }
    3407             :     else
    3408           0 :     { T->t = signe(a) < 0? 3: -1; }
    3409             :   }
    3410           7 :   return 1;
    3411             : }
    3412             : 
    3413             : /* z * sqrt(st_b) / sqrt(st_a) exp(I Pi (t + t0)). Assume that
    3414             :  * sqrt2 = gsqrt(gen_2, prec) or NULL */
    3415             : static GEN
    3416          35 : apply_eta_correction(GEN z, GEN st_a, GEN st_b, GEN t0, GEN sqrt2, long prec)
    3417             : {
    3418          35 :   GEN t, s_a = gel(st_a, 1), s_b = gel(st_b, 1);
    3419             :   cxanalyze_t Ta, Tb;
    3420             :   int ca, cb;
    3421             : 
    3422          35 :   t = gsub(gel(st_b,2), gel(st_a,2));
    3423          35 :   if (t0 != gen_0) t = gadd(t, t0);
    3424          35 :   ca = cxanalyze(&Ta, s_a);
    3425          35 :   cb = cxanalyze(&Tb, s_b);
    3426          35 :   if (ca || cb)
    3427          35 :   { /* compute sqrt(s_b) / sqrt(s_a) in a more efficient way:
    3428             :      * sb = ub sqrt(2)^vb exp(i Pi/4 tb) */
    3429          35 :     GEN u = gdiv(Tb.u,Ta.u);
    3430          35 :     switch(Tb.v - Ta.v)
    3431             :     {
    3432           0 :       case -1: u = gmul2n(u,-1); /* fall through: write 1/sqrt2 = sqrt2/2 */
    3433           7 :       case 1: u = gmul(u, sqrt2? sqrt2: sqrtr_abs(real2n(1, prec)));
    3434             :     }
    3435          35 :     if (!isint1(u)) z = gmul(z, gsqrt(u, prec));
    3436          35 :     t = gadd(t, gmul2n(stoi(Tb.t - Ta.t), -3));
    3437             :   }
    3438             :   else
    3439             :   {
    3440           0 :     z = gmul(z, gsqrt(s_b, prec));
    3441           0 :     z = gdiv(z, gsqrt(s_a, prec));
    3442             :   }
    3443          35 :   return gmul(z, expIPiQ(t, prec));
    3444             : }
    3445             : 
    3446             : /* sqrt(2) eta(2x) / eta(x) */
    3447             : GEN
    3448           7 : weberf2(GEN x, long prec)
    3449             : {
    3450           7 :   pari_sp av = avma;
    3451             :   GEN z, sqrt2, a,b, Ua,Ub, st_a,st_b;
    3452             : 
    3453           7 :   x = upper_to_cx(x, &prec);
    3454           7 :   a = cxredsl2(x, &Ua);
    3455           7 :   b = cxredsl2(gmul2n(x,1), &Ub);
    3456           7 :   if (gequal(a,b)) /* not infrequent */
    3457           0 :     z = gen_1;
    3458             :   else
    3459           7 :     z = gdiv(eta_reduced(b,prec), eta_reduced(a,prec));
    3460           7 :   st_a = eta_correction(a, Ua, 1);
    3461           7 :   st_b = eta_correction(b, Ub, 1);
    3462           7 :   sqrt2 = sqrtr_abs(real2n(1, prec));
    3463           7 :   z = apply_eta_correction(z, st_a, st_b, gen_0, sqrt2, prec);
    3464           7 :   return gerepileupto(av, gmul(z, sqrt2));
    3465             : }
    3466             : 
    3467             : /* eta(x/2) / eta(x) */
    3468             : GEN
    3469          14 : weberf1(GEN x, long prec)
    3470             : {
    3471          14 :   pari_sp av = avma;
    3472             :   GEN z, a,b, Ua,Ub, st_a,st_b;
    3473             : 
    3474          14 :   x = upper_to_cx(x, &prec);
    3475          14 :   a = cxredsl2(x, &Ua);
    3476          14 :   b = cxredsl2(gmul2n(x,-1), &Ub);
    3477          14 :   if (gequal(a,b)) /* not infrequent */
    3478           0 :     z = gen_1;
    3479             :   else
    3480          14 :     z = gdiv(eta_reduced(b,prec), eta_reduced(a,prec));
    3481          14 :   st_a = eta_correction(a, Ua, 1);
    3482          14 :   st_b = eta_correction(b, Ub, 1);
    3483          14 :   z = apply_eta_correction(z, st_a, st_b, gen_0, NULL, prec);
    3484          14 :   return gerepileupto(av, z);
    3485             : }
    3486             : /* exp(-I*Pi/24) * eta((x+1)/2) / eta(x) */
    3487             : GEN
    3488          14 : weberf(GEN x, long prec)
    3489             : {
    3490          14 :   pari_sp av = avma;
    3491             :   GEN z, t0, a,b, Ua,Ub, st_a,st_b;
    3492          14 :   x = upper_to_cx(x, &prec);
    3493          14 :   a = cxredsl2(x, &Ua);
    3494          14 :   b = cxredsl2(gmul2n(gaddgs(x,1),-1), &Ub);
    3495          14 :   if (gequal(a,b)) /* not infrequent */
    3496           7 :     z = gen_1;
    3497             :   else
    3498           7 :     z = gdiv(eta_reduced(b,prec), eta_reduced(a,prec));
    3499          14 :   st_a = eta_correction(a, Ua, 1);
    3500          14 :   st_b = eta_correction(b, Ub, 1);
    3501          14 :   t0 = mkfrac(gen_m1, utoipos(24));
    3502          14 :   z = apply_eta_correction(z, st_a, st_b, t0, NULL, prec);
    3503          14 :   if (typ(z) == t_COMPLEX && isexactzero(real_i(x)))
    3504           0 :     z = gerepilecopy(av, gel(z,1));
    3505             :   else
    3506          14 :     z = gerepileupto(av, z);
    3507          14 :   return z;
    3508             : }
    3509             : GEN
    3510          35 : weber0(GEN x, long flag,long prec)
    3511             : {
    3512          35 :   switch(flag)
    3513             :   {
    3514          14 :     case 0: return weberf(x,prec);
    3515          14 :     case 1: return weberf1(x,prec);
    3516           7 :     case 2: return weberf2(x,prec);
    3517           0 :     default: pari_err_FLAG("weber");
    3518             :   }
    3519             :   return NULL; /* LCOV_EXCL_LINE */
    3520             : }
    3521             : 
    3522             : /* check |q| < 1 */
    3523             : static GEN
    3524          21 : check_unit_disc(const char *fun, GEN q, long prec)
    3525             : {
    3526          21 :   GEN Q = gtofp(q, prec), Qlow;
    3527          21 :   Qlow = (prec > LOWDEFAULTPREC)? gtofp(Q,LOWDEFAULTPREC): Q;
    3528          21 :   if (gcmp(gnorm(Qlow), gen_1) >= 0)
    3529           0 :     pari_err_DOMAIN(fun, "abs(q)", ">=", gen_1, q);
    3530          21 :   return Q;
    3531             : }
    3532             : 
    3533             : GEN
    3534          14 : theta(GEN q, GEN z, long prec)
    3535             : {
    3536             :   long l, n;
    3537          14 :   pari_sp av = avma, av2;
    3538             :   GEN s, c, snz, cnz, s2z, c2z, ps, qn, y, zy, ps2, k, zold;
    3539             : 
    3540          14 :   l = precision(q);
    3541          14 :   n = precision(z); if (n && n < l) l = n;
    3542          14 :   if (l) prec = l;
    3543          14 :   z = gtofp(z, prec);
    3544          14 :   q = check_unit_disc("theta", q, prec);
    3545          14 :   zold = NULL; /* gcc -Wall */
    3546          14 :   zy = imag_i(z);
    3547          14 :   if (gequal0(zy)) k = gen_0;
    3548             :   else
    3549             :   {
    3550           7 :     GEN lq = glog(q,prec); k = roundr(divrr(zy, real_i(lq)));
    3551           7 :     if (signe(k)) { zold = z; z = gadd(z, mulcxmI(gmul(lq,k))); }
    3552             :   }
    3553          14 :   qn = gen_1;
    3554          14 :   ps2 = gsqr(q);
    3555          14 :   ps = gneg_i(ps2);
    3556          14 :   gsincos(z, &s, &c, prec);
    3557          14 :   s2z = gmul2n(gmul(s,c),1); /* sin 2z */
    3558          14 :   c2z = gsubgs(gmul2n(gsqr(c),1), 1); /* cos 2z */
    3559          14 :   snz = s;
    3560          14 :   cnz = c; y = s;
    3561          14 :   av2 = avma;
    3562          14 :   for (n = 3;; n += 2)
    3563         259 :   {
    3564             :     long e;
    3565         273 :     s = gadd(gmul(snz, c2z), gmul(cnz,s2z));
    3566         273 :     qn = gmul(qn,ps);
    3567         273 :     y = gadd(y, gmul(qn, s));
    3568         273 :     e = gexpo(s); if (e < 0) e = 0;
    3569         273 :     if (gexpo(qn) + e < -prec2nbits(prec)) break;
    3570             : 
    3571         259 :     ps = gmul(ps,ps2);
    3572         259 :     c = gsub(gmul(cnz, c2z), gmul(snz,s2z));
    3573         259 :     snz = s; /* sin nz */
    3574         259 :     cnz = c; /* cos nz */
    3575         259 :     if (gc_needed(av,2))
    3576             :     {
    3577           0 :       if (DEBUGMEM>1) pari_warn(warnmem,"theta (n = %ld)", n);
    3578           0 :       gerepileall(av2, 5, &snz, &cnz, &ps, &qn, &y);
    3579             :     }
    3580             :   }
    3581          14 :   if (signe(k))
    3582             :   {
    3583           7 :     y = gmul(y, gmul(powgi(q,sqri(k)),
    3584             :                      gexp(gmul(mulcxI(zold),shifti(k,1)), prec)));
    3585           7 :     if (mod2(k)) y = gneg_i(y);
    3586             :   }
    3587          14 :   return gerepileupto(av, gmul(y, gmul2n(gsqrt(gsqrt(q,prec),prec),1)));
    3588             : }
    3589             : 
    3590             : GEN
    3591           7 : thetanullk(GEN q, long k, long prec)
    3592             : {
    3593             :   long l, n;
    3594           7 :   pari_sp av = avma;
    3595             :   GEN p1, ps, qn, y, ps2;
    3596             : 
    3597           7 :   if (k < 0)
    3598           0 :     pari_err_DOMAIN("thetanullk", "k", "<", gen_0, stoi(k));
    3599           7 :   l = precision(q);
    3600           7 :   if (l) prec = l;
    3601           7 :   q = check_unit_disc("thetanullk", q, prec);
    3602             : 
    3603           7 :   if (!(k&1)) { set_avma(av); return gen_0; }
    3604           7 :   qn = gen_1;
    3605           7 :   ps2 = gsqr(q);
    3606           7 :   ps = gneg_i(ps2);
    3607           7 :   y = gen_1;
    3608           7 :   for (n = 3;; n += 2)
    3609         280 :   {
    3610             :     GEN t;
    3611         287 :     qn = gmul(qn,ps);
    3612         287 :     ps = gmul(ps,ps2);
    3613         287 :     t = gmul(qn, powuu(n, k)); y = gadd(y, t);
    3614         287 :     if (gexpo(t) < -prec2nbits(prec)) break;
    3615             :   }
    3616           7 :   p1 = gmul2n(gsqrt(gsqrt(q,prec),prec),1);
    3617           7 :   if (k&2) y = gneg_i(y);
    3618           7 :   return gerepileupto(av, gmul(p1, y));
    3619             : }
    3620             : 
    3621             : /* q2 = q^2 */
    3622             : static GEN
    3623       13440 : vecthetanullk_loop(GEN q2, long k, long prec)
    3624             : {
    3625       13440 :   GEN ps, qn = gen_1, y = const_vec(k, gen_1);
    3626       13440 :   pari_sp av = avma;
    3627       13440 :   const long bit = prec2nbits(prec);
    3628             :   long i, n;
    3629             : 
    3630       13440 :   if (gexpo(q2) < -2*bit) return y;
    3631       13440 :   ps = gneg_i(q2);
    3632       13440 :   for (n = 3;; n += 2)
    3633       93963 :   {
    3634      107403 :     GEN t = NULL/*-Wall*/, P = utoipos(n), N2 = sqru(n);
    3635      107403 :     qn = gmul(qn,ps);
    3636      107403 :     ps = gmul(ps,q2);
    3637      322209 :     for (i = 1; i <= k; i++)
    3638             :     {
    3639      214806 :       t = gmul(qn, P); gel(y,i) = gadd(gel(y,i), t);
    3640      214806 :       P = mulii(P, N2);
    3641             :     }
    3642      107403 :     if (gexpo(t) < -bit) return y;
    3643       93963 :     if (gc_needed(av,2))
    3644             :     {
    3645           0 :       if (DEBUGMEM>1) pari_warn(warnmem,"vecthetanullk_loop, n = %ld",n);
    3646           0 :       gerepileall(av, 3, &qn, &ps, &y);
    3647             :     }
    3648             :   }
    3649             : }
    3650             : /* [d^i theta/dz^i(q, 0), i = 1, 3, .., 2*k - 1] */
    3651             : GEN
    3652           0 : vecthetanullk(GEN q, long k, long prec)
    3653             : {
    3654           0 :   long i, l = precision(q);
    3655           0 :   pari_sp av = avma;
    3656             :   GEN p1, y;
    3657             : 
    3658           0 :   if (l) prec = l;
    3659           0 :   q = check_unit_disc("vecthetanullk", q, prec);
    3660           0 :   y = vecthetanullk_loop(gsqr(q), k, prec);
    3661           0 :   p1 = gmul2n(gsqrt(gsqrt(q,prec),prec),1);
    3662           0 :   for (i = 2; i <= k; i += 2) gel(y,i) = gneg_i(gel(y,i));
    3663           0 :   return gerepileupto(av, gmul(p1, y));
    3664             : }
    3665             : 
    3666             : /* [d^i theta/dz^i(q, 0), i = 1, 3, .., 2*k - 1], q = exp(2iPi tau) */
    3667             : GEN
    3668           0 : vecthetanullk_tau(GEN tau, long k, long prec)
    3669             : {
    3670           0 :   long i, l = precision(tau);
    3671           0 :   pari_sp av = avma;
    3672             :   GEN q4, y;
    3673             : 
    3674           0 :   if (l) prec = l;
    3675           0 :   if (typ(tau) != t_COMPLEX || gsigne(gel(tau,2)) <= 0)
    3676           0 :     pari_err_DOMAIN("vecthetanullk_tau", "imag(tau)", "<=", gen_0, tau);
    3677           0 :   q4 = expIPiC(gmul2n(tau,-1), prec); /* q^(1/4) */
    3678           0 :   y = vecthetanullk_loop(gpowgs(q4,8), k, prec);
    3679           0 :   for (i = 2; i <= k; i += 2) gel(y,i) = gneg_i(gel(y,i));
    3680           0 :   return gerepileupto(av, gmul(gmul2n(q4,1), y));
    3681             : }
    3682             : 
    3683             : /* Return E_k(tau). Slow if tau is not in standard fundamental domain */
    3684             : GEN
    3685       13727 : cxEk(GEN tau, long k, long prec)
    3686             : {
    3687             :   pari_sp av;
    3688             :   GEN q, y, qn;
    3689       13727 :   long n, b, l = precision(tau);
    3690             : 
    3691       13727 :   if (l) prec = l;
    3692       13727 :   b = bit_accuracy(prec);
    3693             :   /* sum n^(k-1) x^n <= x(1 + (k!-1)x) / (1-x)^k (cf Eulerian polynomials)
    3694             :    * S = \sum_{n > 0} n^(k-1) |q^n/(1-q^n)| <= x(1+(k!-1)x) / (1-x)^(k+1),
    3695             :    * where x = |q| = exp(-2Pi Im(tau)) < 1. Neglegt 2/zeta(1-k) * S if
    3696             :    * (2Pi)^k/(k-1)! x < 2^(-b-1) and k! x < 1. Use log2((2Pi)^k/(k-1)!) < 10 */
    3697       13727 :   if (gcmpgs(imag_i(tau), (M_LN2 / (2*M_PI)) * (b+1+10)) > 0)
    3698          35 :     return real_1(prec);
    3699       13692 :   if (k == 2)
    3700             :   { /* -theta^(3)(tau/2) / theta^(1)(tau/2). Assume that Im tau > 0 */
    3701       13440 :     y = vecthetanullk_loop(qq(tau,prec), 2, prec);
    3702       13440 :     return gdiv(gel(y,2), gel(y,1));
    3703             :   }
    3704         252 :   q = cxtoreal(expIPiC(gneg(gmul2n(tau,1)), prec));
    3705         252 :   av = avma; y = gen_0; qn = q;
    3706         252 :   for(n = 1;; n++)
    3707        2450 :   { /* compute y := sum_{n>0} n^(k-1) / (q^n-1) */
    3708        2702 :     GEN t = gdiv(powuu(n,k-1), gsubgs(qn,1));
    3709        2702 :     if (gequal0(t) || gexpo(t) <= - prec2nbits(prec) - 5) break;
    3710        2450 :     y = gadd(y, t);
    3711        2450 :     qn = gmul(q,qn);
    3712        2450 :     if (gc_needed(av,2))
    3713             :     {
    3714           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"elleisnum");
    3715           0 :       gerepileall(av, 2, &y,&qn);
    3716             :     }
    3717             :   }
    3718         252 :   return gadd(gen_1, gmul(y, gdiv(gen_2, szeta(1-k, prec))));
    3719             : }
    3720             : 
    3721             : /* Lambert W function : solution x of x*exp(x)=y, using Newton. y >= 0 t_REAL.
    3722             :  * Good for low accuracy: the precision remains constant. Not memory-clean */
    3723             : static GEN
    3724         203 : mplambertW0(GEN y)
    3725             : {
    3726         203 :   long bitprec = bit_prec(y) - 2;
    3727             :   GEN x, tmp;
    3728             : 
    3729         203 :   x = mplog(addrs(y,1));
    3730             :   do {
    3731        1176 :    tmp = x;
    3732             :    /* f(x) = log(x)+x-log(y), f'(x) = (1+1/x)
    3733             :     * x *= (1-log(x/y))/(x+1) */
    3734        1176 :     x = mulrr(x, divrr(subsr(1, mplog(divrr(x,y))), addrs(x,1)));
    3735        1176 :   } while (expo(tmp) - expo(subrr(x,tmp)) < bitprec);
    3736         203 :   return x;
    3737             : }
    3738             : 
    3739             : /* Lambert W function using Newton, increasing prec */
    3740             : GEN
    3741         210 : mplambertW(GEN y)
    3742             : {
    3743         210 :   pari_sp av = avma;
    3744             :   GEN x;
    3745         210 :   long p = 1, s = signe(y);
    3746         210 :   ulong mask = quadratic_prec_mask(realprec(y)-1);
    3747             : 
    3748         210 :   if (s<0) pari_err_DOMAIN("Lw", "y", "<", gen_0, y);
    3749         203 :   if(s==0) return rcopy(y);
    3750         203 :   x = mplambertW0(rtor(y, LOWDEFAULTPREC));
    3751         508 :   while (mask!=1)
    3752             :   {
    3753         305 :     p <<= 1; if (mask & 1) p--;
    3754         305 :     mask >>= 1;
    3755         305 :     x = rtor(x, p+2);
    3756         305 :     x = mulrr(x, divrr(subsr(1, mplog(divrr(x,y))),addrs(x,1)));
    3757             :   }
    3758         203 :   return gerepileuptoleaf(av,x);
    3759             : }
    3760             : 
    3761             : /* exp(t (1 + O(t^n))), n >= 1 */
    3762             : static GEN
    3763          91 : serexp0(long v, long n)
    3764             : {
    3765          91 :   GEN y = cgetg(n+3, t_SER), t;
    3766             :   long i;
    3767          91 :   y[1] = evalsigne(1) | evalvarn(v) | evalvalp(0);
    3768          91 :   gel(y,2) = gel(y,3) = gen_1; t = gen_2;
    3769         994 :   for (i = 2; i < n; i++, t = muliu(t,i)) gel(y,i+2) = mkfrac(gen_1,t);
    3770          91 :   gel(y,i+2) = mkfrac(gen_1,t); return y;
    3771             : }
    3772             : 
    3773             : static GEN
    3774          91 : reverse(GEN y)
    3775             : {
    3776          91 :   GEN z = ser2rfrac_i(y);
    3777          91 :   long l = lg(z);
    3778          91 :   return RgX_to_ser(RgXn_reverse(z, l-2), l);
    3779             : }
    3780             : static GEN
    3781         133 : serlambertW(GEN y, long prec)
    3782             : {
    3783             :   GEN x, t, y0;
    3784             :   long n, l, vy, val, v;
    3785             : 
    3786         133 :   if (!signe(y)) return gcopy(y);
    3787         119 :   v = valp(y);
    3788         119 :   vy = varn(y);
    3789         119 :   n = lg(y)-3;
    3790         119 :   y0 = gel(y,2);
    3791         511 :   for (val = 1; val < n; val++)
    3792         483 :     if (!gequal0(polcoef(y, val, vy))) break;
    3793         119 :   if (v < 0) pari_err_DOMAIN("lambertw","valuation", "<", gen_0, y);
    3794         112 :   if (val >= n)
    3795             :   {
    3796          21 :     if (v) return zeroser(vy, n);
    3797          21 :     x = glambertW(y0,prec);
    3798          21 :     return scalarser(x, vy, n+1);
    3799             :   }
    3800          91 :   l = 3 + n/val;
    3801          91 :   if (v)
    3802             :   {
    3803          49 :     t = serexp0(vy, l-3);
    3804          49 :     setvalp(t, 1); /* t exp(t) */
    3805          49 :     t = reverse(t);
    3806             :   }
    3807             :   else
    3808             :   {
    3809          42 :     y = serchop0(y);
    3810          42 :     x = glambertW(y0, prec);
    3811             :     /* (x + t) exp(x + t) = (y0 + t y0/x) * exp(t) */
    3812          42 :     t = gmul(deg1pol_shallow(gdiv(y0,x), y0, vy), serexp0(vy, l-3));
    3813          42 :     t = gadd(x, reverse(serchop0(t)));
    3814             :   }
    3815          91 :   t = gsubst(t, vy, y);
    3816          91 :   return normalize(t);
    3817             : }
    3818             : 
    3819             : GEN
    3820         553 : glambertW(GEN y, long prec)
    3821             : {
    3822             :   pari_sp av;
    3823             :   GEN z;
    3824         553 :   switch(typ(y))
    3825             :   {
    3826         210 :     case t_REAL: return mplambertW(y);
    3827           7 :     case t_COMPLEX: pari_err_IMPL("lambert(t_COMPLEX)");
    3828         336 :     default:
    3829         336 :       av = avma; if (!(z = toser_i(y))) break;
    3830         133 :       return gerepileupto(av, serlambertW(z, prec));
    3831             :   }
    3832         203 :   return trans_eval("lambert",glambertW,y,prec);
    3833             : }
    3834             : 
    3835             : #if 0
    3836             : /* Another Lambert-like function: solution of exp(x)/x=y, y >= e t_REAL,
    3837             :  * using Newton with constant prec. Good for low accuracy */
    3838             : GEN
    3839             : mplambertX(GEN y)
    3840             : {
    3841             :   long bitprec = bit_prec(y)-2;
    3842             :   GEN tmp, x = mplog(y);
    3843             :   if (typ(x) != t_REAL || signe(subrs(x,1))<0)
    3844             :     pari_err(e_MISC,"Lx : argument less than e");
    3845             :   do {
    3846             :     tmp = x;
    3847             :    /* f(x)=x-log(xy), f'(x)=1-1/x */
    3848             :    /* x *= (log(x*y)-1)/(x-1) */
    3849             :     x = mulrr(x, divrr(subrs(mplog(mulrr(x,y)),1), subrs(x,1)));
    3850             :   } while(expo(tmp)-expo(subrr(x,tmp)) < bitprec);
    3851             :   return x;
    3852             : }
    3853             : #endif

Generated by: LCOV version 1.13