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

Generated by: LCOV version 1.11