Code coverage tests

This page documents the degree to which the PARI/GP source code is tested by our public test suite, distributed with the source distribution in directory src/test/. This is measured by the gcov utility; we then process gcov output using the lcov frond-end.

We test a few variants depending on Configure flags on the pari.math.u-bordeaux.fr machine (x86_64 architecture), and agregate them in the final report:

The target is to exceed 90% coverage for all mathematical modules (given that branches depending on DEBUGLEVEL or DEBUGMEM are not covered). This script is run to produce the results below.

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

Generated by: LCOV version 1.13