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.0 lcov report (development 23690-5d6e28857) Lines: 1812 1943 93.3 %
Date: 2019-03-18 05:43:21 Functions: 120 123 97.6 %
Legend: Lines: hit not hit

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

Generated by: LCOV version 1.13