Code coverage tests

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

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

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

LCOV - code coverage report
Current view: top level - basemath - lerch.c (source / functions) Hit Total Coverage
Test: PARI/GP v2.16.2 lcov report (development 29419-8afb0ed749) Lines: 352 365 96.4 %
Date: 2024-07-02 09:03:41 Functions: 23 23 100.0 %
Legend: Lines: hit not hit

          Line data    Source code
       1             : /* Copyright (C) 2022  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; either version 2 of the License, or (at your option) any later
       8             : version. It is distributed in the hope that it will be useful, but WITHOUT
       9             : ANY WARRANTY WHATSOEVER.
      10             : 
      11             : Check the License for details. You should have received a copy of it, along
      12             : with the package; see the file 'COPYING'. If not, write to the Free Software
      13             : Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. */
      14             : 
      15             : #include "pari.h"
      16             : #include "paripriv.h"
      17             : 
      18             : #define DEBUGLEVEL DEBUGLEVEL_trans
      19             : 
      20             : /********************************************************/
      21             : /*                   Hurwitz zeta function              */
      22             : /********************************************************/
      23             : struct Qp_zetahurwitz_t { GEN B, _1, s1; };
      24             : static void
      25         280 : Qp_zetahurwitz_init(struct Qp_zetahurwitz_t *S, long prec, GEN s)
      26             : {
      27         280 :   GEN B, C = gen_1, s1 = gsubgs(s, 1), p = gel(s, 2);
      28         280 :   long j, J = ((equaliu(p,2)? (prec >> 1): prec) + 2) >> 1;
      29         280 :   if (gequal0(s1)) s1 = NULL;
      30         280 :   B = bernvec(J);
      31        2401 :   for (j = 1; j <= J; j++)
      32             :   {
      33        2121 :     GEN t = (j == 1 && !s1)? s: gmul(gaddgs(s, 2*j-3), gaddgs(s, 2*j-2));
      34        2121 :     C = gdivgunextu(gmul(C, t), 2*j-1);
      35        2121 :     gel(B, j+1) = gmul(gel(B, j+1), C); /* B_{2j} * binomial(1-s, 2j) */
      36             :   }
      37         280 :   S->_1 = cvtop(gen_1, p, prec);
      38         280 :   S->s1 = s1;
      39         280 :   S->B = B;
      40         280 : }
      41             : 
      42             : /* v_p(x) < (p==2)?-1: 0; s1 = s-1 or NULL (if s=1) */
      43             : static GEN
      44         952 : Qp_zetahurwitz_0(struct Qp_zetahurwitz_t *S, GEN x)
      45             : {
      46         952 :   GEN z, x2, x2j, s1 = S->s1;
      47         952 :   long j, J = lg(S->B) - 2;
      48             : 
      49         952 :   x = cvtop2(ginv(x), S->_1); z = gmul2n(x, -1);
      50         952 :   z = s1? gmul(s1, z): gadd(Qp_log(x), z);
      51         952 :   x2j = x2 = gsqr(x); z = gaddgs(z, 1);
      52         952 :   for (j = 1;; j++)
      53             :   {
      54        7231 :     z = gadd(z, gmul(gel(S->B, j + 1), x2j));
      55        7231 :     if (j == J) break;
      56        6279 :     x2j = gmul(x2, x2j);
      57             :   }
      58         952 :   if (s1) z = gmul(gdiv(z, s1), Qp_exp(gmul(s1, Qp_log(x))));
      59         952 :   return z;
      60             : }
      61             : /* private (absolute) padicprec */
      62             : static long
      63         518 : pprec(GEN x) { return maxss(valp(x) + precp(x), 1); }
      64             : 
      65             : /* L_p(s, (D, .)); assume s != 1 if D = 1 */
      66             : static GEN
      67          14 : Qp_zeta_i(GEN s, long D)
      68             : {
      69          14 :   pari_sp av = avma;
      70          14 :   GEN z, va, gp = gel(s,2);
      71          14 :   ulong a, p = itou(gp), m;
      72          14 :   long prec = pprec(s);
      73             :   struct Qp_zetahurwitz_t S;
      74             : 
      75          14 :   if (D <= 0) pari_err_DOMAIN("p-adic L-function", "D", "<=", gen_0, stoi(D));
      76          14 :   if (!uposisfundamental(D))
      77           0 :     pari_err_TYPE("p-adic L-function [D not fundamental]", stoi(D));
      78          14 :   Qp_zetahurwitz_init(&S, prec, s);
      79          14 :   m = ulcm(D, p == 2? 4: p); va = coprimes_zv(m);
      80          42 :   for (a = 1, z = gen_0; a <= (m >> 1); a++)
      81          28 :     if (va[a])
      82             :     {
      83          21 :       GEN h = Qp_zetahurwitz_0(&S, uutoQ(a, m));
      84          21 :       if (D != 1 && kross(D, a) < 0) h = gneg(h);
      85          21 :       z = gadd(z, h);
      86             :     }
      87          14 :   z = gdivgs(gmul2n(z, 1), m);
      88          14 :   if (D != 1) z = gmul(z, Qp_exp(gmul(gsubsg(1, s), Qp_log(cvstop2(m, s)))));
      89          14 :   return gerepileupto(av, z);
      90             : }
      91             : GEN
      92          14 : Qp_zeta(GEN s) { return Qp_zeta_i(s, 1); }
      93             : 
      94             : /* s a t_PADIC; gerepileupto-safe. Could be exported */
      95             : static GEN
      96         266 : Qp_zetahurwitz_ii(GEN s, GEN x, long k)
      97             : {
      98         266 :   GEN gp = gel(s,2);
      99         266 :   long p = itou(gp), prec = pprec(s);
     100             :   struct Qp_zetahurwitz_t S;
     101         266 :   Qp_zetahurwitz_init(&S, prec, s);
     102         266 :   if (typ(x) != t_PADIC) x = gadd(x, zeropadic_shallow(gp, prec));
     103         266 :   if (valp(x) >= ((p==2)? -1: 0))
     104             :   {
     105         259 :     GEN z = gen_0;
     106         259 :     long j, M = (p==2)? 4: p;
     107        1561 :     for (j = 0; j < M; j++)
     108             :     {
     109        1302 :       GEN y = gaddsg(j, x);
     110        1302 :       if (valp(y) <= 0)
     111             :       {
     112         924 :         GEN tmp = Qp_zetahurwitz_0(&S, gdivgu(y, M));
     113         924 :         if (k) tmp = gmul(tmp, gpowgs(teich(y), k));
     114         924 :         z = gadd(z, tmp);
     115             :       }
     116             :     }
     117         259 :     return gdivgu(z, M);
     118             :   }
     119           7 :   if (valp(s) < 0) pari_err_DOMAIN("Qp_zetahurwitz", "v(s)", "<", gen_0, s);
     120           7 :   return Qp_zetahurwitz_0(&S, x);
     121             : }
     122             : 
     123             : /* x or s must be p-adic */
     124             : static GEN
     125         266 : Qp_zetahurwitz_i(GEN s, GEN x, long k)
     126             : {
     127         266 :   if (typ(x) == t_PADIC)
     128             :   {
     129         238 :     pari_sp av = avma;
     130         238 :     GEN p = gel(x,2);
     131         238 :     long e = pprec(x);
     132         238 :     e += sdivsi(e, subis(p, 1));
     133         238 :     s = gadd(s, zeropadic_shallow(p, e));
     134         238 :     return gerepileupto(av, Qp_zetahurwitz_ii(s, x, k));
     135             :   }
     136          28 :   return Qp_zetahurwitz_ii(s, x, k);
     137             : }
     138             : 
     139             : GEN
     140         231 : Qp_zetahurwitz(GEN s, GEN x, long k)
     141             : {
     142         231 :   pari_sp av = avma;
     143         231 :   return gerepileupto(av, Qp_zetahurwitz_i(s, x, k));
     144             : }
     145             : 
     146             : static void
     147        8596 : binsplit(GEN *pP, GEN *pR, GEN aN2, GEN isqaN, GEN s, long j, long k, long prec)
     148             : {
     149        8596 :   if (j + 1 == k)
     150             :   {
     151        4347 :     long j2 = j << 1;
     152             :     GEN P;
     153        4347 :     if (!j) P = gdiv(s, aN2);
     154             :     else
     155             :     {
     156        4249 :       P = gmul(gaddgs(s, j2-1), gaddgs(s, j2));
     157        4249 :       P = gdivgunextu(gmul(P, isqaN), j2+1);
     158             :     }
     159        4347 :     if (pP) *pP = P;
     160        4347 :     if (pR) *pR = gmul(bernreal(j2+2, prec), P);
     161             :   }
     162             :   else
     163             :   {
     164             :     GEN P1, R1, P2, R2;
     165        4249 :     binsplit(&P1,pR? &R1: NULL, aN2, isqaN, s, j, (j+k) >> 1, prec);
     166        4249 :     binsplit(pP? &P2: NULL, pR? &R2: NULL, aN2, isqaN, s, (j+k) >> 1, k, prec);
     167        4249 :     if (pP) *pP = gmul(P1,P2);
     168        4249 :     if (pR) *pR = gadd(R1, gmul(P1, R2));
     169             :   }
     170        8596 : }
     171             : 
     172             : /* a0 +  a1 x + O(x^e), e >= 0 */
     173             : static GEN
     174          77 : deg1ser_shallow(GEN a1, GEN a0, long v, long e)
     175          77 : { return RgX_to_ser(deg1pol_shallow(a1, a0, v), e+2); }
     176             : 
     177             : static long
     178         518 : hurwitz_cutoff(GEN s, long bit)
     179             : {
     180         567 :   return typ(s) == t_COMPLEX &&
     181          49 :          fabs(gtodouble(gel(s,2))) > 5.37 * pow(bit, 1.4) / mt_nbthreads();
     182             : }
     183             : 
     184             : /* New zetahurwitz, from Fredrik Johansson. */
     185             : GEN
     186         637 : zetahurwitz(GEN s, GEN x, long der, long bitprec)
     187             : {
     188         637 :   pari_sp av = avma, av2;
     189         637 :   GEN a, ra, ra0, Nx, S1, S2, S3, N2, rx, sch = NULL, s0 = s, x0 = x, y;
     190         637 :   long j, k, m, N, prec0 = nbits2prec(bitprec), prec = prec0, fli = 0;
     191             :   pari_timer T;
     192             : 
     193         637 :   if (der < 0) pari_err_DOMAIN("zetahurwitz", "der", "<", gen_0, stoi(der));
     194         637 :   if (der)
     195             :   {
     196             :     GEN z;
     197          21 :     if (!is_scalar_t(typ(s)))
     198             :     {
     199           7 :       z = deriv(zetahurwitz(s, x, der - 1, bitprec), -1);
     200           7 :       z = gdiv(z, deriv(s, -1));
     201             :     }
     202             :     else
     203             :     {
     204          14 :       if (gequal1(s)) pari_err_DOMAIN("zetahurwitz", "s", "=", gen_1, s0);
     205          14 :       s = deg1ser_shallow(gen_1, s, 0, der+2);
     206          14 :       z = zetahurwitz(s, x, 0, bitprec + der * log2(der));
     207          14 :       z = gmul(mpfact(der), polcoef_i(z, der, -1));
     208             :     }
     209          21 :     return gerepileupto(av,z);
     210             :   }
     211         616 :   if (typ(s) == t_PADIC || typ(x) == t_PADIC)
     212          35 :     return gerepileupto(av, Qp_zetahurwitz_i(s, x, 0));
     213         581 :   switch(typ(x))
     214             :   {
     215         420 :     case t_INT: case t_REAL: case t_FRAC: case t_COMPLEX: break;
     216         161 :     default:
     217         161 :       if (!(y = toser_i(x))) pari_err_TYPE("zetahurwitz", x);
     218         154 :       x = y; x0 = polcoef_i(x, 0, -1); break;
     219             :   }
     220         574 :   rx = grndtoi(real_i(x0), NULL);
     221         574 :   if (typ(rx) != t_INT) pari_err_TYPE("zetahurwitz", x);
     222         574 :   if (x0 == x && signe(rx) <= 0 && gexpo(gsub(x, rx)) < 17 - bitprec)
     223           0 :     pari_err_DOMAIN("zetahurwitz", "x", "<=", gen_0, x);
     224         574 :   switch (typ(s))
     225             :   {
     226             :     long v, pr;
     227         504 :     case t_INT: case t_REAL: case t_FRAC: case t_COMPLEX:
     228         504 :       if (!der && hurwitz_cutoff(s, bitprec))
     229           7 :         return zetahurwitzlarge(s, x, prec);
     230         497 :       break;
     231          70 :     default:
     232          70 :       if (!(y = toser_i(s))) pari_err_TYPE("zetahurwitz", s);
     233          70 :       if (valser(y) < 0) pari_err_DOMAIN("zetahurwitz", "val(s)", "<", gen_0, s);
     234          70 :       s0 = polcoef_i(y, 0, -1);
     235          70 :       switch(typ(s0))
     236             :       {
     237          63 :         case t_INT: case t_REAL: case t_FRAC: case t_COMPLEX: break;
     238           0 :         case t_PADIC: pari_err_IMPL("zetahurwitz(t_SER of t_PADIC)");
     239           7 :         default: pari_err_TYPE("zetahurwitz", s0);
     240             :       }
     241          63 :       sch = gequal0(s0)? y: serchop0(y);
     242          63 :       v = valser(sch);
     243          63 :       pr = (lg(y) + v + 1) / v;
     244          63 :       if (gequal1(s0)) pr += v;
     245          63 :       s = deg1ser_shallow(gen_1, s0, 0, pr);
     246             :     }
     247         560 :   a = gneg(s0); ra = real_i(a); ra0 = ground(ra);
     248         560 :   if (gequal1(s0) && (!sch || gequal0(sch)))
     249          14 :     pari_err_DOMAIN("zetahurwitz", "s", "=", gen_1, s0);
     250         546 :   fli = (gsigne(ra0) >= 0 && gexpo(gsub(a, ra0)) < 17 - bitprec);
     251         546 :   if (!sch && fli)
     252             :   { /* a ~ non negative integer */
     253          14 :     k = itos(gceil(ra)) + 1;
     254          14 :     if (odd(k)) k++;
     255          14 :     N = 1;
     256             :   }
     257             :   else
     258             :   {
     259         532 :     GEN C, ix = imag_i(x0);
     260         532 :     double c = (typ(s) == t_INT)? 1: 20 * log((double)bitprec);
     261         532 :     double rs = gtodouble(ra) + 1;
     262             :     long k0;
     263         532 :     if (fli) a = gadd(a, ghalf); /* hack */
     264         532 :     if (rs > 0)
     265             :     {
     266          49 :       bitprec += (long)ceil(rs * expu(bitprec));
     267          49 :       prec = nbits2prec(bitprec);
     268          49 :       x = gprec_w(x, prec);
     269          49 :       s = gprec_w(s, prec);
     270          49 :       if (sch) sch = gprec_w(sch, prec);
     271             :     }
     272         532 :     k = bitprec * M_LN2 / (1 + dbllambertW0(M_PI / c));
     273         532 :     k0 = itos(gceil(gadd(ra, ghalf))) + 1;
     274         532 :     k = maxss(k0, k);
     275         532 :     if (odd(k)) k++;
     276             :     /* R_k < 2 |binom(a,k+1) B_{k+2}/(k+2)| */
     277         532 :     C = binomial(a, k+1); C = polcoef_i(C, 0, -1);
     278         532 :     C = gmul(C, gdivgu(bernfrac(k+2), k+2));
     279         532 :     C = gmul2n(gabs(C,LOWDEFAULTPREC), bitprec + 1);
     280         532 :     C = gpow(C, ginv(gsubsg(k+1, ra)), LOWDEFAULTPREC);
     281             :     /* need |N + x - 1|^2 > C^2 */
     282         532 :     if (!gequal0(ix))
     283             :     {
     284         133 :       GEN tmp = gsub(gsqr(C), gsqr(ix));
     285         133 :       if (gsigne(tmp) >= 0) C = gsqrt(tmp, LOWDEFAULTPREC);
     286             :     }
     287             :     /* need |N + re(x) - 1| > C */
     288         532 :     C = gceil(gadd(C, gsubsg(1, rx)));
     289         532 :     if (typ(C) != t_INT) pari_err_TYPE("zetahurwitz",s);
     290         532 :     N = signe(C) > 0? itos(C) : 1;
     291         532 :     if (N == 1 && signe(a) > 0)
     292             :     { /* May reduce k if 2Pix > a */
     293             :       /* Need 2 |x^(-K) (B_K/K) binom(a, K-1)| < 2^-bit |x|^-rs |zeta(s,x)|
     294             :        * with K = k+2; N = 1; |zeta(s,x)| ~ |x|^(rs-1);
     295             :        * if a > 0, (B_K/K) binom(a, K-1) < 2 |a / 2Pi|^K */
     296           0 :       double dx = dbllog2(x0), d = 1 + dx + log2(M_PI) - dbllog2(s0);
     297           0 :       if (d > 0)
     298             :       { /* d ~ log2 |2Pi x / a| */
     299           0 :         long K = (long)ceil((bitprec + 1 + dx) / d);
     300           0 :         K = maxss(k0, K);
     301           0 :         if (odd(K)) K++;
     302           0 :         if (K < k) k = K;
     303             :       }
     304             :     }
     305             :   }
     306         546 :   if (gsigne(rx) < 0) N = maxss(N, 1 - itos(rx));
     307         546 :   a = gneg(s);
     308         546 :   if (DEBUGLEVEL>2) timer_start(&T);
     309         546 :   incrprec(prec);
     310         546 :   Nx = gaddsg(N - 1, x);
     311         546 :   Nx = typ(Nx) == t_SER? RgX_gtofp(Nx, prec): gtofp(Nx, prec);
     312         546 :   S1 = S3 = gpow(Nx, a, prec);
     313         546 :   av2 = avma;
     314         546 :   if (gequal1(x)) S1 = dirpowerssum(N, a, 0, prec);
     315             :   else
     316        7532 :     for (m = N - 2; m >= 0; m--)
     317             :     {
     318        7042 :       S1 = gadd(S1, gpow(gaddsg(m,x), a, prec));
     319        7042 :       if ((m & 0xff) == 0) S1 = gerepileupto(av2, S1);
     320             :     }
     321         546 :   if (DEBUGLEVEL>2) timer_printf(&T,"sum from 0 to N - 1");
     322         546 :   constbern(k >> 1);
     323         546 :   N2 = ginv(gsqr(Nx));
     324         546 :   if (typ(s0) == t_INT)
     325             :   {
     326         448 :     S2 = divru(bernreal(k, prec), k);
     327       10479 :     for (j = k - 2; j >= 2; j -= 2)
     328             :     {
     329       10031 :       GEN t = gsubgs(a, j), u = gmul(t, gaddgs(t, 1));
     330       10031 :       u = gmul(gdivgunextu(u, j), gmul(S2, N2));
     331       10031 :       S2 = gadd(divru(bernreal(j, prec), j), u);
     332             :     }
     333         448 :     S2 = gmul(S2, gdiv(a, Nx));
     334             :   }
     335             :   else
     336             :   {
     337          98 :     binsplit(NULL,&S2, gmul2n(Nx,1), N2, s, 0, k >> 1, prec);
     338          98 :     S2 = gneg(S2);
     339             :   }
     340         546 :   S2 = gadd(ghalf, S2);
     341         546 :   if (DEBUGLEVEL>2) timer_printf(&T,"Bernoulli sum");
     342         546 :   S2 = gmul(S3, gadd(gdiv(Nx, gaddsg(1, a)), S2));
     343         546 :   S1 = gprec_wtrunc(gsub(S1, S2), prec0);
     344         546 :   if (sch) return gerepileupto(av, gsubst(S1, 0, sch));
     345         490 :   return gerepilecopy(av, S1);
     346             : }
     347             : 
     348             : /* New Lerch, inspired by Fredrik Johansson. */
     349             : 
     350             : GEN
     351      153717 : lerch_worker(GEN t, GEN E)
     352             : {
     353      153717 :   GEN n, d, T, s = gel(E,1), a = gmul(gel(E,2), t), z = gel(E,3);
     354      153137 :   long p = itos(gel(E,4)), prec = labs(p);
     355      153146 :   d = gadd(gexp(t, prec), z);
     356      153474 :   T = p > 0? t: gneg(t);
     357      153519 :   if (typ(s) == t_INT)
     358       58362 :     n = gmul(gpow(T, s, prec), gexp(a, prec));
     359             :   else /* save one exp */
     360       95157 :     n = gexp(gadd(gmul(s, glog(T, prec)), a), prec);
     361      153754 :   return gdiv(n, d);
     362             : }
     363             : 
     364             : /* tab already computed with N = #tab[1] even */
     365             : static GEN
     366         441 : parintnumgauss(GEN f, GEN a, GEN b, GEN tab, long prec)
     367             : {
     368         441 :   GEN R = gel(tab, 1), W = gel(tab, 2), bma, bpa, S = gen_0, VP, VM, V;
     369         441 :   long n = lg(R) - 1, i, prec2 = prec + EXTRAPREC64;
     370         441 :   a = gprec_wensure(a, prec2);
     371         441 :   b = gprec_wensure(b, prec2);
     372         441 :   VP = cgetg(n + 1, t_VEC); bma = gmul2n(gsub(b, a), -1);
     373         441 :   VM = cgetg(n + 1, t_VEC); bpa = gadd(bma, a);
     374       28203 :   for (i = 1; i <= n; i++)
     375             :   {
     376       27762 :     GEN h = gmul(bma, gel(R, i));
     377       27762 :     gel(VP, i) = gadd(bpa, h);
     378       27762 :     gel(VM, i) = gsub(bpa, h);
     379             :   }
     380         441 :   V = gadd(parapply(f, VP), parapply(f, VM));
     381       28203 :   for (i = 1; i <= n; i++)
     382             :   {
     383       27762 :     S = gadd(S, gmul(gel(W, i), gel(V, i)));
     384       27762 :     S = gprec_wensure(S, prec2);
     385             :   }
     386         441 :   return gprec_wtrunc(gmul(bma, S), prec);
     387             : }
     388             : 
     389             : /* Assume tab computed and a >= 0 */
     390             : static GEN
     391         119 : parintnum(GEN f, GEN a, GEN tab)
     392             : {
     393             :   pari_sp av;
     394         119 :   GEN tabx0 = gel(tab, 2), tabw0 = gel(tab, 3), tabxm = gel(tab, 6);
     395         119 :   GEN tabxp = gel(tab, 4), tabwp = gel(tab, 5), tabwm = gel(tab, 7);
     396         119 :   GEN VP = tabxp, VM = tabxm, x0 = tabx0, S;
     397         119 :   long prec = gprecision(tabw0), L = lg(tabxp), i, fla = 0;
     398         119 :   if (!gequal0(a))
     399             :   {
     400          91 :     if (gexpo(a) <= 0)
     401             :     {
     402          63 :       x0 = gadd(a, x0);
     403       30357 :       for (i = 1; i < L; i++)
     404             :       {
     405       30294 :         gel(VP, i) = gadd(a, gel(VP, i));
     406       30294 :         gel(VM, i) = gadd(a, gel(VM, i));
     407             :       }
     408             :     }
     409             :     else
     410             :     {
     411          28 :       x0 = gmul(a, gaddsg(1, x0)); fla = 1;
     412        5404 :       for (i = 1; i < L; i++)
     413             :       {
     414        5376 :         gel(VP, i) = gmul(a, gaddsg(1, gel(VP, i)));
     415        5376 :         gel(VM, i) = gmul(a, gaddsg(1, gel(VM, i)));
     416             :       }
     417             :     }
     418             :   }
     419         119 :   VP = parapply(f, VP);
     420         119 :   VM = parapply(f, VM); av = avma;
     421         119 :   S = gmul(tabw0, closure_callgen1(f, x0));
     422       49257 :   for (i = 1; i < L; i++)
     423             :   {
     424       49138 :     S = gadd(S, gadd(gmul(gel(tabwp, i), gel(VP, i)),
     425       49138 :                      gmul(gel(tabwm, i), gel(VM, i))));
     426       49138 :     if ((i & 0x7f) == 1) S = gerepileupto(av, S);
     427       49138 :     S = gprec_wensure(S, prec);
     428             :   }
     429         119 :   if (fla) S = gmul(S, a);
     430         119 :   return gmul(S, gel(tab, 1));
     431             : }
     432             : 
     433             : static GEN
     434          84 : refine(GEN A)
     435             : {
     436          84 :   long n = lg(A) - 1, i;
     437          84 :   GEN B = cgetg(2 * n, t_VEC);
     438         231 :   for (i = 1; i < n; i++)
     439             :   {
     440         147 :     gel(B, 2 * i - 1) = gel(A, i);
     441         147 :     gel(B, 2 * i) = gmul2n(gadd(gel(A, i), gel(A, i + 1)), -1);
     442             :   }
     443          84 :   gel(B, 2 * n - 1) = gel(A, n); return B;
     444             : }
     445             : 
     446             : /* Here L = [a1, a2, a3,...] integration vertices. Refine by splitting
     447             :  * intervals. */
     448             : static GEN
     449          84 : parintnumgaussadapt(GEN f, GEN L, GEN tab, long bit)
     450             : {
     451          84 :   GEN Rold = gen_0, Rnew;
     452          84 :   long i, ct = 0, prec = nbits2prec(bit);
     453         168 :   while (ct <= 5)
     454             :   {
     455         168 :     Rnew = gen_0;
     456         609 :     for (i = 1; i < lg(L) - 1; i++)
     457         441 :       Rnew = gadd(Rnew, parintnumgauss(f, gel(L, i), gel(L, i + 1), tab, prec));
     458         168 :     if (ct && gexpo(gsub(Rnew, Rold)) - gexpo(Rnew) < 10 - bit) return Rnew;
     459          84 :     ct++; Rold = Rnew; L = refine(L);
     460             :   }
     461           0 :   if (DEBUGLEVEL) err_printf("intnumgaussadapt: possible accuracy loss");
     462           0 :   return Rnew; /* Possible accuracy loss */
     463             : }
     464             : 
     465             : /* Here b = [oo, r], so refine by increasing integration step m */
     466             : static GEN
     467          42 : parintnumadapt(GEN f, GEN a, GEN b, GEN tab, long bit)
     468             : {
     469          42 :   GEN Rold = gen_0, Rnew;
     470          42 :   long m = 0, prec = nbits2prec(bit);
     471          42 :   if (!tab) tab = intnuminit(gen_0, b, 0, prec);
     472         119 :   while (m <= 5)
     473             :   {
     474         119 :     Rnew = parintnum(f, a, tab);
     475         119 :     if (m && gexpo(gsub(Rnew, Rold)) - gexpo(Rnew) < 10 - bit) return Rnew;
     476          77 :     m++; Rold = Rnew; tab = intnuminit(gen_0, b, m, prec);
     477             :   }
     478           0 :   if (DEBUGLEVEL) err_printf("intnumadapt: possible accuracy loss");
     479           0 :   return Rnew; /* Possible accuracy loss */
     480             : }
     481             : 
     482             : static int
     483         210 : iscplx(GEN z) { long t = typ(z); return is_real_t(t) || t == t_COMPLEX; }
     484             : 
     485             : static GEN
     486          14 : lerch_easy(GEN z, GEN s, GEN a, long B)
     487             : {
     488          14 :   long n, prec = nbits2prec(B + 32);
     489          14 :   GEN zn, ms = gneg(s), S = gpow(a, ms, prec);
     490          14 :   zn = z = gtofp(z, prec);
     491        3808 :   for (n = 1;; n++, zn = gmul(zn, z))
     492             :   {
     493        3808 :     S = gadd(S, gmul(zn, gpow(gaddgs(a, n), ms, prec)));
     494        3808 :     if (gexpo(zn) <= - B - 5) return S;
     495             :   }
     496             : }
     497             : 
     498             : static GEN
     499         112 : _lerchphi(GEN z, GEN s, GEN a, long prec)
     500             : {
     501         112 :   GEN res = NULL, L, LT, J, rs, mleft, left, right, top, w, Linf, tabg;
     502             :   GEN E, f, fm;
     503         112 :   long B = prec2nbits(prec), MB = 3 - B, NB, prec2;
     504             :   entree *ep;
     505             : 
     506         112 :   if (gexpo(z) < MB) return gpow(a, gneg(s), prec);
     507         112 :   if (gexpo(gsubgs(z, 1)) < MB) return zetahurwitz(s, a, 0, B); /* z ~ 1 */
     508         112 :   if (gexpo(gaddgs(z, 1)) < MB) /* z ~ -1 */
     509             :   {
     510           7 :     GEN tmp = gsub(zetahurwitz(s, gmul2n(a, -1), 0, B),
     511             :                    zetahurwitz(s, gmul2n(gaddgs(a, 1), -1), 0, B));
     512           7 :     return gmul(gpow(gen_2, gneg(s), prec), tmp);
     513             :   }
     514         105 :   if (gcmpgs(gmulsg(10, gabs(z, prec)), 9) <= 0) /* |z| <= 9/10 */
     515          14 :     return lerch_easy(z, s, a, B);
     516          91 :   if (gcmpgs(real_i(a), 2) < 0)
     517          49 :     return gadd(gpow(a, gneg(s), prec),
     518             :                 gmul(z, _lerchphi(z, s, gaddgs(a, 1), prec)));
     519          42 :   NB = (long)ceil(B + M_PI * fabs(gtodouble(imag_i(s))));
     520          42 :   prec2 = nbits2prec(NB);
     521          42 :   z = gprec_w(z, prec2); /* |z| > 9/10 */
     522          42 :   s = gprec_w(s, prec2);
     523          42 :   a = gprec_w(a, prec2); /* Re(a) >= 2 */
     524          42 :   rs = ground(real_i(s)); L = glog(z, prec2); /* Re(L) > -0.11 */
     525          42 :   ep = is_entry("_lerch_worker");
     526          42 :   E = mkvec4(gsubgs(s, 1), gsubsg(1, a), gneg(z), stoi(prec2));
     527          42 :   f = snm_closure(ep, mkvec(E));
     528          42 :   E = shallowcopy(E); gel(E,4) = stoi(-prec2);
     529          42 :   fm = snm_closure(ep, mkvec(E));
     530          42 :   Linf = mkvec2(mkoo(), real_i(a));
     531          42 :   if (gexpo(gsub(s, rs)) < MB && gcmpgs(rs, 1) >= 0)
     532             :   { /* s ~ positive integer */
     533          14 :     if (gcmp(gabs(imag_i(L), prec2), sstoQ(1, 4)) < 0 && gsigne(real_i(L)) >= 0)
     534           7 :     { /* Re(L) >= 0, |Im(L)| < 1/4 */
     535           7 :       GEN t = gsigne(imag_i(z)) > 0 ? gen_m1: gen_1;
     536           7 :       GEN LT1 = gaddgs(gabs(L, prec2), 1);
     537           7 :       LT = mkvec4(gen_0, mkcomplex(gen_0, t), mkcomplex(LT1, t), LT1);
     538           7 :       tabg = intnumgaussinit(2*(NB >> 2) + 60, prec2);
     539           7 :       J = parintnumgaussadapt(f, LT, tabg, NB);
     540           7 :       J = gadd(J, parintnumadapt(f, LT1, Linf, NULL, NB));
     541             :     }
     542           7 :     else J = parintnumadapt(f, gen_0, Linf, NULL, NB);
     543          14 :     return gdiv(J, ggamma(s, prec2));
     544             :   }
     545          28 :   tabg = intnumgaussinit(2*(NB >> 2) + 60, prec2);
     546          28 :   if (gcmp(gabs(imag_i(L), prec2), ghalf) > 0) /* |Im(L)| > 1/2 */
     547          14 :     left = right = top = gmin(gmul2n(gabs(imag_i(L), prec2), -1), gen_1);
     548             :   else
     549             :   {
     550          14 :     res = gdiv(gpow(gneg(L), s, prec2), gmul(L, gpow(z, a, prec2)));
     551          14 :     left = gaddgs(gmax(gen_0, gneg(real_i(L))), 1);
     552          14 :     top = gaddgs(gabs(imag_i(L), prec2), 1);
     553          14 :     right = gaddgs(gabs(L, prec2), 1);
     554             :   }
     555          28 :   w = expIPiC(gsubgs(s, 1), prec2);
     556          28 :   mleft = gneg(left);
     557          28 :   if (gexpo(imag_i(z)) < MB && gexpo(imag_i(a)) < MB && gexpo(imag_i(s)) < MB
     558           7 :       && gcmpgs(real_i(z), 1) < 0)
     559             :   { /* (z, s, a) real, z < 1 */
     560           7 :     LT = mkvec3(right, mkcomplex(right, top), mkcomplex(mleft, top));
     561           7 :     J = imag_i(gdiv(parintnumgaussadapt(f, LT, tabg, NB), w));
     562           7 :     LT = mkvec2(mkcomplex(mleft, top), mleft);
     563           7 :     J = gmul2n(gadd(J, imag_i(parintnumgaussadapt(fm, LT, tabg, NB))), 1);
     564           7 :     J = mulcxI(J);
     565             :   }
     566             :   else
     567             :   {
     568          21 :     GEN mtop = gneg(top);
     569          21 :     LT = mkvec3(right, mkcomplex(right, top), mkcomplex(mleft, top));
     570          21 :     J = gdiv(parintnumgaussadapt(f, LT, tabg, NB), w);
     571          21 :     LT = mkvec2(mkcomplex(mleft, top), mkcomplex(mleft, mtop));
     572          21 :     J = gadd(J, parintnumgaussadapt(fm, LT, tabg, NB));
     573          21 :     LT = mkvec3(mkcomplex(mleft, mtop), mkcomplex(right, mtop), right);
     574          21 :     J = gadd(J, gmul(parintnumgaussadapt(f, LT, tabg, NB), w));
     575             :   }
     576          28 :   J = gadd(J, gmul(gsub(w, ginv(w)), parintnumadapt(f, right, Linf, NULL, NB)));
     577          28 :   J = gdiv(J, PiI2(prec2)); if (res) J = gadd(J, res);
     578          28 :   return gneg(gmul(ggamma(gsubsg(1, s), prec2), J));
     579             : }
     580             : /* lerchphi(z,-k,a)=
     581             :  *  -1/(z-1)*sum(q=0,k,(z/(z-1))^q*sum(j=0,q,(-1)^j*(j+a)^k*binomial(q,j)))
     582             :  * zetahurwitz(-k,a)=-B(k+1,a)/(k+1) */
     583             : GEN
     584          56 : lerchphi(GEN z, GEN s, GEN a, long prec)
     585             : {
     586          56 :   pari_sp av = avma;
     587          56 :   if (!iscplx(z)) pari_err_TYPE("lerchphi", z);
     588          56 :   if (!iscplx(s)) pari_err_TYPE("lerchphi", s);
     589          56 :   if (!iscplx(a)) pari_err_TYPE("lerchphi", a);
     590          56 :   return gerepileupto(av, _lerchphi(z, s, a, prec));
     591             : }
     592             : 
     593             : GEN
     594          14 : lerchzeta(GEN s, GEN a, GEN lam, long prec)
     595             : {
     596          14 :   pari_sp av = avma;
     597          14 :   GEN z = gexp(gmul(PiI2(prec), lam), prec);
     598          14 :   if (!iscplx(z)) pari_err_TYPE("lerchzeta", z);
     599          14 :   if (!iscplx(s)) pari_err_TYPE("lerchzeta", s);
     600          14 :   if (!iscplx(a)) pari_err_TYPE("lerchzeta", a);
     601          14 :   if (hurwitz_cutoff(s, prec)) return lerchzetalarge(s, a, lam, prec);
     602           7 :   return gerepileupto(av, _lerchphi(z, s, a, prec));
     603             : }

Generated by: LCOV version 1.16