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-bordeaux1.fr machine (x86_64 architecture), and agregate them in the final report:

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

LCOV - code coverage report
Current view: top level - basemath - trans3.c (source / functions) Hit Total Coverage
Test: PARI/GP v2.8.0 lcov report (development 16358-a11f489) Lines: 1669 1827 91.4 %
Date: 2014-04-11 Functions: 106 111 95.5 %
Legend: Lines: hit not hit | Branches: + taken - not taken # not executed Branches: 821 1105 74.3 %

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

Generated by: LCOV version 1.9