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.13.0 lcov report (development 25825-edc109b529) Lines: 1905 2037 93.5 %
Date: 2020-09-21 06:08:33 Functions: 131 134 97.8 %
Legend: Lines: hit not hit

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

Generated by: LCOV version 1.13