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 - lfunlarge.c (source / functions) Hit Total Coverage
Test: PARI/GP v2.18.0 lcov report (development 29589-e347ff5c25) Lines: 479 518 92.5 %
Date: 2024-10-06 09:07:08 Functions: 73 75 97.3 %
Legend: Lines: hit not hit

          Line data    Source code
       1             : /* Copyright (C) 2000  The PARI group.
       2             : 
       3             : This file is part of the PARI/GP package.
       4             : 
       5             : PARI/GP is free software; you can redistribute it and/or modify it under the
       6             : terms of the GNU General Public License as published by the Free Software
       7             : Foundation. It is distributed in the hope that it will be useful, but WITHOUT
       8             : ANY WARRANTY WHATSOEVER.
       9             : 
      10             : Check the License for details. You should have received a copy of it, along
      11             : with the package; see the file 'COPYING'. If not, write to the Free Software
      12             : Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. */
      13             : 
      14             : #include "pari.h"
      15             : #include "paripriv.h"
      16             : 
      17             : /*****************************************************************/
      18             : /*             Program to compute L(chi,s)                       */
      19             : /*      for Im(s) large, chi primitive Dirichlet character       */
      20             : /*      In the present branch, only Tyagi's method is used      */
      21             : /*****************************************************************/
      22             : /*
      23             :    In addition, C can also be a polynomial defining an abelian
      24             :    extension of Q.
      25             : */
      26             : 
      27             : /*****************************************************************/
      28             : /*                      Character programs                       */
      29             : /*****************************************************************/
      30             : /* A character, always assumed primitive can be given in the following formats:
      31             :  * - omitted or 0: special to zetaRS,
      32             :  * - a t_INT: assumed to be a discriminant,
      33             :  * - a t_INTMOD: a conrey character,
      34             :  * - a pair [G,chi] or [bnr,chi],
      35             :  * - [C1,C2,...]~ where the Ci are characters as above with same moduli. */
      36             : 
      37             : /* Given a list of linit/ldata for chars of same conductor F, return
      38             :  * [Vecan, F, Parities, Gaussums] */
      39             : static GEN
      40        1085 : mycharinit(GEN C, long bit)
      41             : {
      42             :   GEN L, LVC, LE, LGA;
      43        1085 :   long F = 0, i, j, lc = lg(C), prec;
      44             : 
      45        1085 :   bit += 64; prec = nbits2prec(bit);
      46        1085 :   L = cgetg(lc, t_VEC);
      47        1085 :   LE = cgetg(lc, t_VECSMALL);
      48        1085 :   LGA = cgetg(lc, t_VEC);
      49        2205 :   for (i = 1; i < lc; i++)
      50             :   {
      51        1120 :     GEN gv, ga, gm, ro, ldata = gel(C, i);
      52             :     long e;
      53        1120 :     if (is_linit(ldata)) ldata = linit_get_ldata(ldata);
      54        1120 :     gv = ldata_get_gammavec(ldata); e = itou(gel(gv, 1));
      55        1120 :     gm = ldata_get_conductor(ldata);
      56        1120 :     ro = ldata_get_rootno(ldata);
      57        1120 :     if (isintzero(ro)) ro = lfunrootno(ldata, bit);
      58        1120 :     ga = gmul(ro, gsqrt(gm, prec)); if (e) ga = mulcxI(ga);
      59        1120 :     gel(LGA, i) = ga;
      60        1120 :     LE[i] = e;
      61        1120 :     if (i == 1) F = itos(gm); /* constant */
      62        1120 :     gel(L, i) = lfunan(ldata, F, prec);
      63             :   }
      64        1085 :   if (lc == 2 && is_vec_t(typ(gmael(L,1,1))))
      65             :   { /* multichar */
      66           7 :     LGA = gel(LGA,1); lc = lg(LGA);
      67           7 :     LVC = gel(L,1);
      68           7 :     LE = const_vecsmall(lc-1, LE[1]); /* FIXME: can handle mixed values */
      69             :   }
      70             :   else
      71             :   {
      72        1078 :     LVC = cgetg(F + 1, t_VEC);
      73        4459 :     for (j = 1; j <= F; j++)
      74             :     {
      75        3381 :       GEN v = cgetg(lc, t_VEC);
      76        7049 :       for (i = 1; i < lc; i++) gel(v, i) = gmael(L, i, j);
      77        3381 :       gel(LVC, j) = v;
      78             :     }
      79             :   }
      80        1085 :   return mkvec5(LVC, stoi(F), LE, LGA, grootsof1(2*F, prec));
      81             : }
      82             : 
      83             : /* n >= 1 and #VC = F, the conductor of the character or multicharacter X.
      84             :  * VC contains [X(1),X(2),...X(F)] */
      85             : static GEN
      86      618472 : mycall(GEN VC, long n)
      87             : {
      88      618472 :   long F = lg(VC) - 1;
      89      618472 :   GEN z = n <= F ? gel(VC, n) : gel(VC, ((n - 1) % F) + 1);
      90      618472 :   return gequal0(z)? NULL: z;
      91             : }
      92             : 
      93      187267 : static GEN get_chivec(GEN VCALL) { return gel(VCALL, 1); }
      94      763801 : static long get_modulus(GEN VCALL) { return itos(gel(VCALL, 2)); }
      95      185965 : static GEN get_signat(GEN VCALL) { return gel(VCALL, 3); }
      96          21 : static GEN get_gauss(GEN VCALL) { return gel(VCALL, 4); }
      97      184663 : static GEN get_chiZ(GEN VCALL) { return gel(VCALL, 5); }
      98             : 
      99             : /* (-1)^A[i] * conj(B[i]) */
     100             : static GEN
     101           7 : gnegconj(GEN A, GEN B)
     102             : {
     103           7 :   long i, l = lg(A);
     104           7 :   GEN W = cgetg(l, t_VEC);
     105          14 :   for (i = 1; i < l; i++)
     106           7 :   { GEN b = gconj(gel(B,i)); gel(W,i) = A[i]? gneg(b): b; }
     107           7 :   return W;
     108             : }
     109             : /* g(conj(CHI)) */
     110             : static GEN
     111           7 : gaussconj(GEN VCALL)
     112           7 : { return gnegconj(get_signat(VCALL), get_gauss(VCALL)); }
     113             : 
     114             : static GEN
     115           7 : myinitconj(GEN VCALL)
     116             : {
     117           7 :   GEN CONJ = shallowcopy(VCALL);
     118           7 :   gel(CONJ, 1) = gconj(get_chivec(VCALL));
     119           7 :   gel(CONJ, 4) = gaussconj(VCALL); return CONJ;
     120             : }
     121             : 
     122             : /********************************************************************/
     123             : /*                          Driver Program                          */
     124             : /********************************************************************/
     125             : 
     126             : /* assume |Im(s)| >> 1, in particular s is not a negative integer */
     127             : static GEN
     128          14 : applyfuneq(GEN gau, GEN s, GEN z, long odd, long q, long bitprec)
     129             : {
     130             :   GEN t, S;
     131             :   long prec;
     132          14 :   if (!gequal0(s)) bitprec += maxss(gexpo(s), 0);
     133          14 :   prec = nbits2prec(bitprec);
     134          14 :   if (odd) gau = mulcxmI(gau);
     135          14 :   S = gmul(Pi2n(-1, prec), gsubgs(s, odd));
     136          14 :   t = ginv(gmul2n(gmul(gcos(S, prec), ggamma(s, prec)), 1));
     137          14 :   t = gmul(gpow(gdivgs(Pi2n(1, prec), q), s, prec), t);
     138          14 :   return gmul(gmul(gau, t), z);
     139             : }
     140             : 
     141             : static GEN RZchi(GEN VCALL, GEN s, long prec);
     142             : 
     143             : /* VCALL already initialized */
     144             : static GEN
     145        1085 : lfunlarge_char(GEN VCALL, GEN s, long bitprec)
     146             : {
     147        1085 :   pari_sp av = avma;
     148             :   GEN sig, tau, z;
     149        1085 :   long funeq = 0, ts = typ(s), stau, flconj, q;
     150        1085 :   if (!is_real_t(ts) && ts != t_COMPLEX) pari_err_TYPE("lfunlarge_char", s);
     151        1085 :   sig = real_i(s); tau = imag_i(s);
     152        1085 :   if (gexpo(tau) < 1) pari_err_DOMAIN("lfun","im(s)", "<", gen_2, tau);
     153        1085 :   stau = gsigne(tau);
     154        1085 :   if (stau < 0) { tau = gneg(tau); VCALL = myinitconj(VCALL); }
     155        1085 :   if (gcmp(sig, ghalf) < 0) { funeq = 1; sig = gsubsg(1, sig); }
     156        1085 :   flconj = ((stau > 0 && funeq) || (stau < 0 && !funeq));
     157        1085 :   q = get_modulus(VCALL); bitprec += gexpo(stoi(q));
     158        1085 :   z = RZchi(VCALL, mkcomplex(sig, tau), nbits2prec(bitprec));
     159        1085 :   if (flconj) z = gconj(z);
     160        1085 :   if (funeq)
     161             :   {
     162          14 :     GEN odd = get_signat(VCALL), gau = get_gauss(VCALL), Vz;
     163          14 :     long lC = lg(gau), j;
     164          14 :     Vz = cgetg(lC, t_VEC);
     165          28 :     for (j = 1; j < lC; j++)
     166          14 :       gel(Vz,j) = applyfuneq(gel(gau,j), s, gel(z,j), odd[j], q, bitprec);
     167          14 :     z = Vz;
     168             :   }
     169        1085 :   return gerepilecopy(av, z);
     170             : }
     171             : 
     172             : static GEN
     173          14 : lfungetchars(GEN pol)
     174             : {
     175          14 :   GEN L, F, v, bnf = Buchall(pol_x(1), 0, LOWDEFAULTPREC);
     176             :   GEN w, condall, bnr;
     177             :   long i, l, lc;
     178          14 :   condall = rnfconductor(bnf, pol); bnr = gel(condall, 2);
     179          14 :   L = bnrchar(bnr, gel(condall, 3), NULL); lc = lg(L);
     180          14 :   F = cgetg(lc, t_VEC);
     181          77 :   for (i = 1; i < lc; i++)
     182             :   {
     183          63 :     GEN chi = gel(L, i), cond = bnrconductor_raw(bnr, chi);
     184          63 :     gel(F, i) = gcoeff(gel(cond,1), 1, 1);
     185             :   }
     186          14 :   w = vec_equiv(F); l = lg(w); v = cgetg(l, t_COL);
     187          42 :   for (i = 1; i < l; i++)
     188             :   {
     189          28 :     GEN wi = gel(w, i), vi;
     190          28 :     long j, li = lg(wi);
     191          28 :     gel(v,i) = vi = cgetg(li, t_VEC);
     192          28 :     if (li == 2 && equali1(gel(F, wi[1]))) /* common conductor is 1 */
     193          14 :       gel(vi,1) = lfunmisc_to_ldata_shallow(gen_1);
     194             :     else
     195             :     {
     196          63 :       for (j = 1; j < li; j++)
     197          49 :         gel(vi,j) = lfunmisc_to_ldata_shallow(mkvec2(bnr, gel(L, wi[j])));
     198             :     }
     199             :   }
     200          14 :   return v;
     201             : }
     202             : 
     203             : /********************************************************************/
     204             : /*            NEW RS IMPLEMENTATION FROM SANDEEP TYAGI              */
     205             : /********************************************************************/
     206             : /* See arXiv:2203.02509v2 */
     207             : 
     208        5523 : static long m_n0(GEN sel) { return itos(gel(sel, 1)); }
     209      576072 : static GEN m_r0(GEN sel) { return gel(sel, 2); }
     210        1085 : static GEN m_al(GEN sel) { return gel(sel, 3); }
     211      576072 : static GEN m_aleps(GEN sel) { return gel(sel, 4); }
     212        4396 : static GEN m_h(GEN sel) { return gel(sel, 5); }
     213        1113 : static GEN m_lin(GEN sel) { return gel(sel, 6); }
     214        1113 : static long m_np(GEN sel) { return itou(gel(sel, 7)); }
     215             : 
     216             : static GEN
     217        3297 : phi_hat(GEN x, long prec)
     218             : {
     219             :   GEN y;
     220        3297 :   if (gsigne(imag_i(x)) > 0)
     221        1078 :     y = gneg(gexpm1(gneg(gmul(PiI2(prec), x)), prec));
     222             :   else
     223        2219 :     y = gexpm1(gmul(PiI2(prec), x), prec);
     224        3297 :   return ginv(y);
     225             : }
     226             : 
     227             : static GEN
     228        3241 : phi_hat_h0(GEN sel, long k, long prec)
     229             : {
     230        3241 :   GEN t = gdiv(gsubsg(m_n0(sel) + k, m_r0(sel)), m_aleps(sel));
     231        3241 :   return phi_hat(gdiv(gasinh(t, prec), m_h(sel)), prec);
     232             : }
     233             : 
     234             : /* v[i] = A[i] * (a + (-1)^E[i] b) */
     235             : static GEN
     236      612774 : mul_addsub(GEN A, GEN a, GEN b, GEN E)
     237             : {
     238      612774 :   long i, l = lg(E);
     239      612774 :   GEN v = cgetg(l, t_VEC);
     240     1308996 :   for (i = 1; i < l; i++)
     241      696222 :     gel(v,i) = gmul(gel(A,i), E[i]? gsub(a, b): gadd(a, b));
     242      612774 :   return v;
     243             : }
     244             : 
     245             : static GEN
     246      184663 : wd(GEN VCALL, GEN pmd, GEN x, long prec)
     247             : {
     248      184663 :   GEN VC = get_chivec(VCALL), E = get_signat(VCALL), Z = get_chiZ(VCALL);
     249      184663 :   GEN ex, emx, xpmd = gmul(x, pmd), y = NULL;
     250      184663 :   long md = get_modulus(VCALL), N = 2*md, k;
     251      184663 :   ex = gexp(mulcxI(xpmd), prec); emx = ginv(ex);
     252      797437 :   for (k = 1; k <= (md-1) / 2; k++)
     253             :   {
     254      612774 :     GEN xc = mycall(VC, k);
     255      612774 :     if (xc)
     256             :     {
     257      612774 :       GEN p3, p2, p1 = gmul(xc, gel(Z, Fl_neg(Fl_sqr(k,N), N) + 1));
     258      612774 :       GEN a = gmul(ex, gel(Z, N - k + 1)), b = gmul(emx, gel(Z, k + 1));
     259      612774 :       GEN c = gmul(ex, gel(Z, k + 1)), d = gmul(emx, gel(Z, N - k + 1));
     260      612774 :       if (odd(md))
     261             :       {
     262      594831 :         p2 = ginv(mulcxmI(gmul2n(gsub(a,b), -1))); /* 1 / sin(xpmd - kpmd) */
     263      594831 :         p3 = ginv(mulcxmI(gmul2n(gsub(c,d), -1))); /* 1 / sin(xpmd + kpmd) */
     264             :       }
     265             :       else
     266             :       {
     267       17943 :         p2 = mulcxI(gdiv(gadd(a,b), gsub(a,b))); /* cotan(xpmd - kpmd) */
     268       17943 :         p3 = mulcxI(gdiv(gadd(c,d), gsub(c,d))); /* cotan(xpmd + kpmd) */
     269             :       }
     270      612774 :       p1 = mul_addsub(p1, p2, p3, E);
     271      612774 :       y = y ? gadd(y, p1) : p1;
     272             :     }
     273             :   }
     274      184663 :   return mulcxmI(gdivgs(y, N));
     275             : }
     276             : 
     277             : static GEN
     278        1085 : series_h0(long n0, GEN s, GEN VCALL, long fl, long prec)
     279             : {
     280        1085 :   GEN C = get_modulus(VCALL) == 1? NULL: get_chivec(VCALL);
     281        1085 :   GEN R = pardirpowerssumfun(C, n0, gneg(s), fl, prec);
     282        1085 :   if (C) return R;
     283         686 :   if (fl) return mkvec2(mkvec(gel(R,1)), mkvec(gel(R,2)));
     284         658 :   return mkvec(R);
     285             : }
     286             : 
     287             : static GEN
     288        1113 : series_residues_h0(GEN sel, GEN s, GEN VCALL, long prec)
     289             : {
     290        1113 :   long n0 = m_n0(sel), np = m_np(sel), k;
     291        1113 :   GEN val = gen_0, VC = get_chivec(VCALL);
     292        4424 :   for (k = maxss(1 - np, 1 - n0); k <= 1 + np; k++)
     293             :   {
     294        3311 :     GEN nk = mycall(VC, n0 + k); /* n0 + k > 0 */
     295        3311 :     if (nk) val = gadd(val, gmul(gmul(phi_hat_h0(sel, k, prec), nk),
     296             :                                  gpow(stoi(n0 + k), gneg(s), prec)));
     297             :   }
     298        1113 :   return val;
     299             : }
     300             : 
     301             : static GEN
     302      572831 : integrand_h0(GEN sel, GEN s, GEN VCALL, GEN x, long prec)
     303             : {
     304      572831 :   pari_sp av = avma;
     305      572831 :   long md = get_modulus(VCALL);
     306      572831 :   GEN r0 = m_r0(sel), aleps = m_aleps(sel), zn, p1;
     307      572831 :   GEN pmd = divru(mppi(prec), md), ix = ginv(x);
     308      572831 :   zn = gadd(r0, gdivgs(gmul(aleps, gsub(x, ix)), 2));
     309      572831 :   p1 = gmul(expIxy(pmd, gsqr(zn), prec),
     310             :             gmul(gpow(zn, gneg(s), prec), gmul(aleps, gadd(x, ix))));
     311      572831 :   if (md == 1)
     312      388168 :     p1 = gdiv(mkvec(mulcxI(p1)), gmul2n(gsin(gmul(pmd, zn), prec), 2));
     313             :   else
     314      184663 :     p1 = gdivgs(gmul(p1, wd(VCALL, pmd, zn, prec)), -2);
     315      572831 :   return gerepileupto(av, p1);
     316             : }
     317             : 
     318             : static GEN
     319        1113 : integral_h0(GEN sel, GEN s, GEN VCALL, long prec)
     320             : {
     321        1113 :   GEN lin_grid = m_lin(sel), S = gen_0;
     322        1113 :   pari_sp av = avma;
     323        1113 :   long j, l = lg(lin_grid);
     324      569604 :   for (j = 1; j < l; j++)
     325             :   {
     326      568491 :     S = gadd(S, integrand_h0(sel, s, VCALL, gel(lin_grid, j), prec));
     327      568491 :     if ((j & 0xff) == 0) S = gerepileupto(av, S);
     328             :   }
     329        1113 :   return gerepileupto(av, gmul(m_h(sel), S));
     330             : }
     331             : 
     332             : /* a + log |x|, a t_REAL, low accuracy */
     333             : static GEN
     334        4508 : addrlogabs(GEN a, GEN x)
     335             : {
     336        4508 :   long prec = DEFAULTPREC;
     337        4508 :   if (gequal0(x)) return real_m2n(64, prec); /* -oo */
     338        4438 :   switch(typ(x))
     339             :   {
     340        4438 :     case t_COMPLEX: return gmul2n(glog(cxnorm(x), prec), -1);
     341           0 :     case t_REAL: break;
     342           0 :     default: x = gtofp(x, prec);
     343             :   }
     344           0 :   return addrr(a, logr_abs(x));
     345             : }
     346             : 
     347             : struct fun_q_t { GEN sel, s, VCALL, B; };
     348             : static GEN
     349        4340 : fun_q(void *E, GEN x)
     350             : {
     351        4340 :   struct fun_q_t *S = (struct fun_q_t *)E;
     352        4340 :   GEN z = integrand_h0(S->sel,S->s,S->VCALL, gexp(x,DEFAULTPREC), DEFAULTPREC);
     353        4340 :   if (typ(z) == t_VEC) z = vecsum(z);
     354        4340 :   return addrlogabs(S->B, z);
     355             : }
     356             : static GEN
     357        2254 : brent_q(void *E, GEN (*f)(void*,GEN), GEN q_low, GEN q_hi)
     358             : {
     359        2254 :   if (gsigne(f(E, q_low)) * gsigne(f(E, q_hi)) >= 0) return NULL;
     360           0 :   return zbrent(E, f, q_low, q_hi, LOWDEFAULTPREC);
     361             : }
     362             : static GEN
     363        1127 : findq(void *E, GEN (*f)(void*,GEN), GEN lq, long B)
     364             : {
     365        1127 :   GEN q_low, q_hi, q_right, q_left, q_est = gasinh(lq, LOWDEFAULTPREC);
     366        1127 :   q_low = gdivgs(gmulsg(4, q_est), 5);
     367        1127 :   q_hi = gdivgs(gmulsg(3, q_est), 2);
     368        1127 :   q_right = brent_q(E, f, q_low, q_hi); if (!q_right) q_right = q_est;
     369        1127 :   q_left = brent_q(E, f, gneg(q_low), gneg(q_hi)); if (!q_left) q_left = q_est;
     370        1127 :   return bitprecision0(gmax(q_right, q_left), B);
     371             : }
     372             : 
     373             : static GEN
     374        1085 : set_q_value(GEN sel, GEN s, GEN VCALL, long prec)
     375             : {
     376             :   struct fun_q_t E;
     377        1085 :   GEN al = m_al(sel), lq;
     378        1085 :   long md = get_modulus(VCALL), LD = DEFAULTPREC;
     379        1085 :   E.sel = sel; E.s = s; E.VCALL = VCALL, E.B = mulur(prec, mplog2(LD));
     380        1085 :   lq = gdiv(gsqrt(gdiv(gmulsg(md, E.B), Pi2n(1, LD)), LD), al);
     381        1085 :   return findq((void*)&E, &fun_q, lq, prec);
     382             : }
     383             : 
     384             : static GEN
     385        1099 : setlin_grid_exp(GEN h, long m, long prec)
     386             : {
     387        1099 :   GEN w, vex = gpowers(gexp(h, prec), (m - 1)/2);
     388             :   long i;
     389        1099 :   w = cgetg(m+1, t_VEC); gel(w, (m + 1)/2) = gen_1;
     390      283108 :   for (i = (m + 3)/2; i <= m; i++)
     391             :   {
     392      282009 :     GEN t1 = gel(vex, i - ((m - 1)/2));
     393      282009 :     gel(w, i) = t1; gel(w, (m + 1) - i) = ginv(t1);
     394             :   }
     395        1099 :   return w;
     396             : }
     397             : 
     398             : static long
     399        1099 : get_m(GEN q, long prec)
     400             : {
     401        1099 :   GEN t = divrr(mulur(4 * prec2nbits(prec), mplog2(prec)), sqrr(mppi(prec)));
     402        1099 :   return 2*itos(gfloor(mulrr(q, t))) + 1;
     403             : }
     404             : 
     405             : static GEN
     406        1085 : RZinit(GEN s, GEN VCALL, GEN numpoles, long prec)
     407             : {
     408             :   GEN sel, al, aleps, n0, r0, q, h;
     409        1085 :   long md = get_modulus(VCALL), m;
     410        1085 :   al = gcmpgs(gabs(imag_i(s), prec), 100) < 0 ? ginv(stoi(4)) : gen_1;
     411        1085 :   r0 = gsqrt(gdiv(gmulgs(s, md), PiI2(prec)), prec);
     412        1085 :   n0 = gfloor(gsub(real_i(r0), imag_i(r0)));
     413        1085 :   aleps = gmul(al, gexp(PiI2n(-2, prec), prec));
     414        1085 :   sel = mkvecn(7, n0, r0, al, aleps, NULL, NULL, numpoles);
     415        1085 :   q = set_q_value(sel, s, VCALL, prec);
     416        1085 :   m = get_m(q, prec);
     417        1085 :   gel(sel,5) = h = divru(q, (m - 1) >> 1);
     418        1085 :   gel(sel,6) = setlin_grid_exp(h, m, prec);
     419        1085 :   return sel;
     420             : }
     421             : 
     422             : static GEN
     423         441 : xpquo_one(GEN s, GEN cs, GEN ga, long odd, long md, long prec)
     424             : {
     425         441 :   GEN rho, a = odd? gen_1: gen_0, z = divsr(md, mppi(prec));
     426         441 :   rho = gmul(gdiv(gpow(gen_I(), gdivgs(gneg(a), 2), prec), gsqrt(ga, prec)),
     427             :              gpow(stoi(md), ginv(stoi(4)), prec));
     428         441 :   return gmul(gdiv(gconj(gmul(rho, gpow(z, gdivgs(cs, 2), prec))),
     429             :                    gmul(rho, gpow(z, gdivgs(s, 2), prec))),
     430             :               gexp(gsub(gconj(glngamma(gdivgs(gadd(cs, a), 2), prec)),
     431             :                         glngamma(gdivgs(gadd(s, a), 2), prec)), prec));
     432             : }
     433             : 
     434             : static GEN
     435        1085 : xpquo(GEN s, GEN VCALL, long prec)
     436             : {
     437        1085 :   long md = get_modulus(VCALL), n, lve, i;
     438        1085 :   GEN cd = NULL, z, pz, cs, VC = get_chivec(VCALL), ve, R;
     439        1085 :   if (!gequal0(s)) prec = nbits2prec(prec2nbits(prec) + maxss(gexpo(s), 0));
     440        1085 :   z = gexp(gdivgs(PiI2(prec), -md), prec);
     441        1085 :   if (md == 1)
     442         686 :     return gmul(gpow(mppi(prec), gsub(s, ghalf), prec),
     443             :                 gexp(gsub(glngamma(gdivgs(gsubsg(1, s), 2), prec),
     444             :                           glngamma(gdivgs(s, 2), prec)), prec));
     445         399 :   pz = gpowers(z, md - 1);
     446        2786 :   for (n = 1; n < md; n++)
     447             :   {
     448        2387 :     GEN xn = mycall(VC, n);
     449        2387 :     if (xn)
     450             :     {
     451        2352 :       GEN tmp = gmul(xn, gel(pz, n + 1));
     452        2352 :       cd = cd ? gadd(cd, tmp) : tmp;
     453             :     }
     454             :   }
     455         399 :   cs = gsubsg(1, gconj(s));
     456         399 :   ve = get_signat(VCALL); lve = lg(ve); R = cgetg(lve, t_VEC);
     457         840 :   for (i = 1; i < lve; i++)
     458         441 :     gel(R, i) = xpquo_one(s, cs, gel(cd, i), ve[i], md, prec);
     459         399 :   if (lve == 2) R = gel(R, 1);
     460         399 :   return R;
     461             : }
     462             : 
     463             : static GEN
     464        1113 : total_value(GEN serh0, GEN sel, GEN s, GEN VCALL, long prec)
     465             : {
     466        1113 :   return gadd(integral_h0(sel, s, VCALL, prec),
     467             :               gsub(serh0, series_residues_h0(sel, s, VCALL, prec)));
     468             : }
     469             : static GEN
     470        1085 : dirichlet_ours(GEN s, GEN VCALL, long prec)
     471             : {
     472        1085 :   int fl = !gequal(real_i(s), ghalf);
     473        1085 :   GEN sel = RZinit(s, VCALL, gen_1, prec);
     474        1085 :   GEN S1, S2, serh0 = series_h0(m_n0(sel), s, VCALL, fl, prec);
     475        1085 :   if (!fl)
     476        1057 :     S2 = S1 = total_value(serh0, sel, s, VCALL, prec);
     477             :   else
     478             :   {
     479          28 :     S1 = total_value(gel(serh0,1), sel, s, VCALL, prec);
     480          28 :     S2 = total_value(gel(serh0,2), sel, gsubsg(1, gconj(s)), VCALL, prec);
     481             :   }
     482        1085 :   return gadd(S1, vecmul(xpquo(s, VCALL, prec), gconj(S2)));
     483             : }
     484             : 
     485             : /* assume |Im(s)| > 2^-bitprec */
     486             : static GEN
     487        1085 : RZchi(GEN VCALL, GEN s, long prec)
     488             : {
     489        1085 :   long prec2 = prec + EXTRAPREC64;
     490        1085 :   return gprec_wtrunc(dirichlet_ours(gprec_w(s, prec2), VCALL, prec2), prec);
     491             : }
     492             : 
     493             : /********************************************************************/
     494             : /*                         Utility Functions                        */
     495             : /********************************************************************/
     496             : /* lam = 0, return L(s); else Lambda(s) */
     497             : static GEN
     498        1085 : lfuncharall(GEN VCALL, GEN s, long lam, long bitprec)
     499             : {
     500        1085 :   GEN ve, P, Q, R, z = lfunlarge_char(VCALL, s, bitprec);
     501             :   long l, i, q, prec;
     502        1085 :   if (!lam) return z;
     503         882 :   ve = get_signat(VCALL); l = lg(ve);
     504         882 :   q = get_modulus(VCALL); prec = nbits2prec(bitprec);
     505         882 :   R = cgetg(l, t_VEC);
     506         882 :   Q = divur(q, mppi(prec));
     507         882 :   P = (q == 1 || zv_equal0(ve))? NULL: gsqrt(utoipos(q), prec);
     508        1764 :   for (i = 1; i < l; i++)
     509             :   {
     510         882 :     GEN se = gmul2n(gaddgs(s, ve[i]), -1), r;
     511         882 :     if (lam == 1)
     512             :     {
     513           0 :       r = gmul(gpow(Q, se, prec), ggamma(se, prec));
     514           0 :       if (P && ve[i]) r = gdiv(r, P);
     515             :     }
     516             :     else
     517             :     {
     518         882 :       r = gadd(gmul(se, glog(Q, prec)), glngamma(se, prec));
     519         882 :       if (P && ve[i]) r = gsub(r, glog(P, prec));
     520             :     }
     521         882 :     gel(R, i) = r;
     522             :   }
     523         882 :   return lam == 1 ? vecmul(R, z) : gadd(R, glog(z, prec));
     524             : }
     525             : 
     526             : static GEN
     527          28 : lfunlargeall_from_chars(GEN v, GEN s, long lam, long bit)
     528             : {
     529          28 :   long i, l = lg(v);
     530          84 :   for (i = 1; i < l; i++)
     531             :   {
     532          56 :     GEN w = mycharinit(gel(v, i), bit), L = lfuncharall(w, s, lam, bit);
     533          56 :     gel(v, i) = lam==-1 ? vecsum(L): vecprod(L);
     534             :   }
     535          28 :   return lam==-1 ? vecsum(v): vecprod(v);
     536             : }
     537             : static GEN
     538        1057 : lfunlargeall(GEN ldata, GEN s, long lam, long bit)
     539             : {
     540             :   GEN w, an;
     541        1057 :   if (lg(ldata) == 2)
     542             :   { /* HACK: ldata[1] a t_DESC_PRODUCT from lfunabelianrelinit / Q */
     543          14 :     GEN v = lfunprod_get_fact(linit_get_tech(gel(ldata,1)));
     544             :     long i, l;
     545          14 :     v = shallowcopy(gel(v,1)); l = lg(v);
     546          42 :     for (i = 1; i < l; i++) gel(v,i) = mkvec(gel(v,i));
     547          14 :     return lfunlargeall_from_chars(v, s, lam, bit);
     548             :   }
     549        1043 :   an = gel(ldata_get_an(ldata), 2);
     550        1043 :   switch(ldata_get_type(ldata))
     551             :   {
     552          14 :     case t_LFUN_NF:
     553             :     {
     554          14 :       GEN v = lfungetchars(nf_get_pol(an));
     555          14 :       return lfunlargeall_from_chars(v, s, lam, bit);
     556             :     }
     557           7 :     case t_LFUN_CHIGEN:
     558             :     {
     559           7 :       GEN chi = gmael(an, 2, 2);
     560           7 :       if (lg(chi) > 1 && is_vec_t(typ(gel(chi,1))))
     561             :       { /* multi char */
     562           7 :         w = mycharinit(mkcol(ldata), bit);
     563           7 :         return lfuncharall(w, s, lam, bit);
     564             :       }
     565             :     }
     566             :     default: /* single char */
     567        1022 :       w = mycharinit(mkcol(ldata), bit);
     568        1022 :       return gel(lfuncharall(w, s, lam, bit), 1);
     569             :   }
     570             : }
     571             : 
     572             : GEN
     573         189 : lfunlarge(GEN CHI, GEN s, long bit)
     574         189 : { return lfunlargeall(CHI, s, 0, bit); }
     575             : 
     576             : GEN
     577           0 : lfunlambdalarge(GEN CHI, GEN s, long bit)
     578           0 : { return lfunlargeall(CHI, s, 1, bit); }
     579             : 
     580             : GEN
     581         868 : lfunloglambdalarge(GEN CHI, GEN s, long bit)
     582         868 : { return lfunlargeall(CHI, s, -1, bit); }
     583             : 
     584             : /********************************************************************/
     585             : /*           LERCH RS IMPLEMENTATION FROM SANDEEP TYAGI             */
     586             : /********************************************************************/
     587             : 
     588             : static GEN
     589          56 : double_exp_residue_pos_h(GEN selsm, long k, long ind, long prec)
     590             : {
     591          56 :   long nk = itos(gel(selsm, 1)) + k;
     592          56 :   GEN r = gel(selsm, 2), ale = gel(selsm, 3), aor = gel(selsm, 4);
     593          56 :   GEN h = gel(selsm, 5), t = gen_0;
     594          56 :   switch(ind)
     595             :   {
     596          28 :     case 0: t = gaddsg(nk, aor); break;
     597           0 :     case 1: t = gneg(gaddsg(nk, aor)); break;
     598          28 :     case 2: t = gsubsg(nk, aor); break;
     599             :   }
     600          56 :   return gdiv(gasinh(gdiv(gsub(t, r), ale), prec), h);
     601             : }
     602             : 
     603             : static GEN
     604          56 : phi_hat_h(GEN selsm, long m, long ind, long prec)
     605          56 : { return phi_hat(double_exp_residue_pos_h(selsm, m, ind, prec), prec); }
     606             : 
     607             : static long
     608          56 : myex(GEN is)
     609          56 : { return gequal0(is) ? 0 : maxss(0, 2 + gexpo(is)); }
     610             : 
     611             : static GEN
     612          56 : gaminus(GEN s, long prec)
     613             : {
     614          56 :   GEN is = imag_i(s), tmp;
     615          56 :   long B = prec2nbits(prec);
     616          56 :   if (gcmpgs(is, -5*B) < 0) return gen_0;
     617          56 :   prec = nbits2prec(B + myex(is));
     618          56 :   tmp = gexp(gsub(glngamma(s, prec), gmul(PiI2n(-1, prec), s)), prec);
     619          56 :   return bitprecision0(tmp, B);
     620             : }
     621             : 
     622             : static GEN
     623          28 : gaplus(GEN s, long prec)
     624             : {
     625          28 :   GEN is = imag_i(s), tmp;
     626          28 :   long B = prec2nbits(prec);
     627          28 :   if (gcmpgs(is, 5*B) > 0) return gen_0;
     628           0 :   prec = nbits2prec(B + myex(is));
     629           0 :   tmp = gexp(gadd(glngamma(s, prec), gmul(PiI2n(-1, prec), s)), prec);
     630           0 :   return bitprecision0(tmp, B);
     631             : }
     632             : 
     633             : GEN
     634        3514 : serh_worker(GEN k, GEN z, GEN a, GEN ns, GEN gprec)
     635             : {
     636        3514 :   long prec = itou(gprec);
     637        3514 :   return gmul(gpow(z, k, prec), gpow(gadd(a, k), ns, prec));
     638             : }
     639             : 
     640             : static void
     641          28 : set_arg(GEN worker, GEN z, GEN a, GEN ns, long prec)
     642          28 : { gel(worker, 7) = mkvec4(z, a, ns, stoi(prec)); }
     643             : 
     644             : 
     645             : static GEN
     646          14 : series_h0l(GEN worker, long n0, GEN s, GEN a, GEN lam, long prec)
     647             : {
     648          14 :   GEN z = typ(lam) == t_INT ? gen_1 : gexp(gmul(PiI2(prec), lam), prec);
     649          14 :   set_arg(worker, z, a, gneg(s), prec);
     650          14 :   return parsum(gen_0, utoi(n0), worker);
     651             : }
     652             : 
     653             : static GEN
     654          14 : series_h1(GEN worker, long n1, GEN s, GEN a, GEN lam, long prec)
     655             : {
     656          14 :   GEN mP, pre_factor, z, sn = gsubgs(s, 1);
     657          14 :   GEN ini = gequal0(lam) ? gen_1 : gen_0;
     658          14 :   pre_factor = gaplus(gneg(sn), prec);
     659          14 :   if (gequal0(pre_factor)) return gen_0;
     660           0 :   mP = gneg(PiI2(prec));
     661           0 :   pre_factor = gmul(gmul(pre_factor, gexp(gmul(mP, gmul(a, lam)), prec)),
     662             :                     gpow(Pi2n(1, prec), sn, prec));
     663           0 :   z = typ(a) == t_INT ? gen_1 : gexp(gmul(mP, a), prec);
     664           0 :   set_arg(worker, z, lam, sn, prec);
     665           0 :   return gmul(pre_factor,  parsum(ini, stoi(n1 - 1), worker));
     666             : }
     667             : 
     668             : static GEN
     669          14 : series_h2(GEN worker, long n2, GEN s, GEN a, GEN lam, long prec)
     670             : {
     671          14 :   GEN P, pre_factor, z, sn = gsubgs(s, 1);
     672          14 :   pre_factor = gaminus(gneg(sn), prec);
     673          14 :   if (gequal0(pre_factor)) return gen_0;
     674          14 :   P = PiI2(prec);
     675          14 :   pre_factor = gmul(gmul(pre_factor, gexp(gmul(gneg(P), gmul(a, lam)), prec)),
     676             :                     gpow(Pi2n(1, prec), sn, prec));
     677          14 :   z = typ(a) == t_INT ? gen_1 : gexp(gmul(P, a), prec);
     678          14 :   set_arg(worker, z, gneg(lam), sn, prec);
     679          14 :   return gmul(pre_factor, parsum(gen_1, stoi(n2), worker));
     680             : }
     681             : 
     682             : static GEN
     683          14 : series_residues_h0l(long numpoles, GEN selsm0, GEN s, GEN a, GEN lam, long prec)
     684             : {
     685          14 :   GEN P, val = gen_0, ra = real_i(a);
     686          14 :   long n0 = m_n0(selsm0), k;
     687          14 :   P = PiI2(prec);
     688          42 :   for (k = -numpoles + 1; k <= numpoles; k++)
     689          28 :     if (gsigne(gaddsg(n0 + k, ra)) > 0)
     690          28 :       val = gadd(val, gmul(gmul(phi_hat_h(selsm0, k, 0, prec),
     691             :                                 gexp(gmul(P, gmulgs(lam, n0 + k)), prec)),
     692             :                            gpow(gaddsg(n0 + k, a), gneg(s), prec)));
     693          14 :   return val;
     694             : }
     695             : 
     696             : static GEN
     697          14 : series_residues_h1(long numpoles, GEN selsm1, GEN s, GEN a, GEN lam, long prec)
     698             : {
     699          14 :   GEN mP, val = gen_0, rlam = real_i(lam), pre_factor, sn = gsubgs(s, 1);
     700          14 :   long n1 = m_n0(selsm1), k;
     701          14 :   pre_factor = gaplus(gneg(sn), prec);
     702          14 :   if (gequal0(pre_factor)) return gen_0;
     703           0 :   mP = gneg(PiI2(prec));
     704           0 :   pre_factor = gmul(gmul(pre_factor, gexp(gmul(mP, gmul(a, lam)), prec)),
     705             :                     gpow(Pi2n(1, prec), sn, prec));
     706           0 :   for (k = -numpoles; k <= numpoles - 1; k++)
     707           0 :     if (gsigne(gaddsg(n1 + k, rlam)) > 0)
     708           0 :       val = gadd(val, gmul(gmul(phi_hat_h(selsm1, k, 1, prec),
     709             :                                 gexp(gmul(mP, gmulgs(a, n1 + k)), prec)),
     710             :                            gpow(gaddsg(n1 + k, lam), sn, prec)));
     711           0 :   return gmul(pre_factor, val);
     712             : }
     713             : 
     714             : static GEN
     715          14 : series_residues_h2(long numpoles, GEN selsm2, GEN s, GEN a, GEN lam, long prec)
     716             : {
     717          14 :   GEN P, val = gen_0, rlam = real_i(lam), pre_factor, sn = gsubgs(s, 1);
     718          14 :   long n2 = m_n0(selsm2), k;
     719          14 :   pre_factor = gaminus(gneg(sn), prec);
     720          14 :   if (gequal0(pre_factor)) return gen_0;
     721          14 :   P = PiI2(prec);
     722          14 :   pre_factor = gmul(gmul(pre_factor, gexp(gmul(gneg(P), gmul(a, lam)), prec)),
     723             :                     gpow(Pi2n(1, prec), sn, prec));
     724          42 :   for (k = -numpoles + 1; k <= numpoles; k++)
     725          28 :     if (gsigne(gsubsg(n2 + k, rlam)) > 0)
     726          28 :       val = gsub(val, gmul(gmul(phi_hat_h(selsm2, k, 2, prec),
     727             :                                 gexp(gmul(P, gmulgs(a, n2 + k)), prec)),
     728             :                            gpow(gsubsg(n2 + k, lam), sn, prec)));
     729          14 :   return gmul(pre_factor, val);
     730             : }
     731             : 
     732             : static GEN
     733        3430 : integrand_h0l(GEN selsm0, GEN s, GEN alam1, GEN x, long prec)
     734             : {
     735        3430 :   GEN r0 = gel(selsm0, 2), ale = gel(selsm0, 3), a = gel(selsm0, 4);
     736        3430 :   GEN ix = ginv(x), zn = gadd(r0, gmul2n(gmul(ale, gsub(x, ix)), -1));
     737        3430 :   GEN P = PiI2n(0, prec), den, num;
     738        3430 :   den = gexpm1(gmul(P, gmul2n(gsub(zn,a), 1)), prec);
     739        3430 :   num = gexp(gmul(gmul(P, zn), gsub(alam1, zn)), prec);
     740        3430 :   num = gmul(gmul(gmul(num, ale), gmul2n(gadd(x, ix), -1)),
     741             :              gpow(zn, gneg(s), prec));
     742        3430 :   return gdiv(num, den);
     743             : }
     744             : 
     745             : static GEN
     746        6860 : integrand_h12(GEN selsm1, GEN s, GEN alam1, GEN x, long prec)
     747             : {
     748        6860 :   GEN r1 = gel(selsm1, 2), ale = gel(selsm1, 3), lam = gel(selsm1, 4);
     749        6860 :   GEN ix = ginv(x), zn = gadd(r1, gmul2n(gmul(ale, gsub(x, ix)), -1));
     750        6860 :   GEN P = PiI2n(0, prec), den, num, y;
     751        6860 :   den = gexpm1(gmul(P, gmul2n(gadd(zn,lam), 1)), prec);
     752        6860 :   num = gexp(gmul(gmul(P, zn), gadd(alam1, zn)), prec);
     753        6860 :   num = gmul(gmul(gmul(num, ale), gmul2n(gadd(x, ix), -1)),
     754             :              gpow(zn, gsubgs(s, 1), prec));
     755        6860 :   y = gdiv(num, den);
     756        6860 :   if (gcmp(garg(zn, prec), Pi2n(-2, prec)) > 0)
     757        1722 :     y = gmul(y, gexp(gmul(PiI2(prec), gsubsg(1, s)), prec));
     758        6860 :   return y;
     759             : }
     760             : 
     761             : static GEN
     762          14 : integral_h0l(GEN lin_grid, GEN selsm0, GEN s, GEN a, GEN lam, long prec)
     763             : {
     764          14 :   GEN A = gaddgs(gmul2n(gadd(a, lam),1), 1), S = gen_0;
     765          14 :   pari_sp av = avma;
     766          14 :   long j, l = lg(lin_grid);
     767             : 
     768        3388 :   for (j = 1; j < l; j++)
     769             :   {
     770        3374 :     S = gadd(S, integrand_h0l(selsm0, s, A, gel(lin_grid, j), prec));
     771        3374 :     if ((j & 0xff) == 0) S = gerepileupto(av, S);
     772             :   }
     773          14 :   S = gmul(m_h(selsm0), S);
     774          14 :   A = gmul(a, gaddsg(1, gadd(a, gmul2n(lam, 1))));
     775          14 :   return gmul(S, gexp(gneg(gmul(PiI2n(0, prec), A)), prec));
     776             : }
     777             : 
     778             : /* do not forget a minus sign for index 2 */
     779             : static GEN
     780          28 : integral_h12(GEN lin_grid, GEN selsm1, GEN s, GEN a, GEN lam, long prec)
     781             : {
     782          28 :   GEN A, E, S = gen_0, ga = gaminus(gsubsg(1, s), prec);
     783          28 :   pari_sp av = avma;
     784          28 :   long j, l = lg(lin_grid);
     785             : 
     786          28 :   if (gequal0(ga)) return S;
     787          28 :   A = gaddgs(gmul2n(gadd(a,lam), 1), 1);
     788        6776 :   for (j = 1; j < l; j++)
     789             :   {
     790        6748 :     S = gadd(S, integrand_h12(selsm1, s, A, gel(lin_grid, j), prec));
     791        6748 :     if ((j & 0xff) == 0) S = gerepileupto(av, S);
     792             :   }
     793          28 :   if (gequal0(S)) return gen_0;
     794          28 :   S = gmul(m_h(selsm1), S);
     795          28 :   E = gexp(gmul(PiI2n(0, prec), gmul(lam, gaddgs(lam, 1))), prec);
     796          28 :   return gmul(gmul(gmul(S, ga), E), gpow(Pi2n(1, prec), gsubgs(s, 1), prec));
     797             : }
     798             : 
     799             : struct _fun_q0_t { GEN sel, s, alam1, B; };
     800             : static GEN
     801          56 : _fun_q0(void *E, GEN x)
     802             : {
     803          56 :   struct _fun_q0_t *S = (struct _fun_q0_t*)E;
     804          56 :   GEN z = integrand_h0l(S->sel, S->s, S->alam1, x, DEFAULTPREC);
     805          56 :   return addrlogabs(S->B, z);
     806             : }
     807             : static GEN
     808         112 : _fun_q12(void *E, GEN x)
     809             : {
     810         112 :   struct _fun_q0_t *S = (struct _fun_q0_t*)E;
     811         112 :   GEN z = integrand_h12(S->sel, S->s, S->alam1, x, DEFAULTPREC);
     812         112 :   return addrlogabs(S->B, z);
     813             : }
     814             : 
     815             : static GEN
     816          14 : RZLERinit(GEN s, GEN a, GEN lam, GEN al, GEN numpoles, long prec)
     817             : {
     818             :   GEN eps, r0, r1, r2, h, lin_grid, q, q0, q1, q2, sel0, sel1, sel2, lq;
     819          14 :   GEN pinv = ginv(PiI2(prec)), c = gmul2n(gadd(a, lam), -1), n0, n1, n2, c2;
     820             :   long m;
     821             :   struct _fun_q0_t E;
     822             : 
     823          14 :   if (!al || gequal0(al))
     824           0 :     al = gcmpgs(gabs(imag_i(s), prec), 100) < 0 ? ginv(stoi(4)) : gen_1;
     825          14 :   c2 = gsub(gsqr(c), gmul(s, pinv));
     826          14 :   r0 = gadd(c, gsqrt(c2, prec));
     827          14 :   r1 = gsqrt(gadd(c2, pinv), prec);
     828          14 :   r2 = gsub(r1, c);
     829          14 :   r1 = gneg(gadd(r1, c));
     830          14 :   n0 = gfloor(gsub(gadd(real_i(r0), imag_i(r0)), a));
     831          14 :   n1 = gneg(gfloor(gadd(gsub(real_i(r1), imag_i(r1)), real_i(lam))));
     832          14 :   n2 = gfloor(gadd(gsub(real_i(r2), imag_i(r2)), real_i(lam)));
     833             : 
     834          14 :   E.s = s; E.alam1 = gaddgs(gmul2n(gadd(a, lam), 1), 1);
     835          14 :   E.B = mulur(prec, mplog2(prec));
     836          14 :   lq = gmul(al, sqrtr_abs(mulrr(divsr(prec, Pi2n(1, DEFAULTPREC)),
     837             :                                 mplog2(DEFAULTPREC))));
     838          14 :   eps = gexp(PiI2n(-2, prec), prec);
     839          14 :   E.sel = sel0 = mkvec5(n0, r0, gdiv(al, eps), a, gen_0);
     840          14 :   q0 = findq(&E, &_fun_q0, lq, prec);
     841             : 
     842          14 :   if (!gequal1(al)) lq = gdiv(lq, gsqr(al));
     843          14 :   E.sel = sel1 = mkvec5(n1, r1, gmul(al, eps), lam, gen_0);
     844          14 :   q1 = findq(&E, &_fun_q12, lq, prec);
     845          14 :   E.sel = sel2 = mkvec5(n2, r2, gmul(al, eps), lam, gen_0);
     846          14 :   q2 = findq(&E, &_fun_q12, lq, prec);
     847          14 :   q = vecmax(mkvec3(q0, q1, q2)); m = get_m(q, prec);
     848          14 :   gel(sel0, 5) = gel(sel1, 5) = gel(sel2, 5) = h = divru(q, (m-1) >> 1);
     849          14 :   lin_grid = setlin_grid_exp(h, m, prec);
     850          14 :   if (!numpoles) numpoles = gen_1;
     851          14 :   return mkvec5(sel0, sel1, sel2, lin_grid, numpoles);
     852             : }
     853             : 
     854          42 : static GEN add3(GEN x, GEN y, GEN z) { return gadd(x, gadd(y,z)); }
     855          14 : static GEN addsub(GEN x, GEN y, GEN z) { return gadd(x, gsub(y,z)); }
     856             : 
     857             : static GEN
     858          14 : lerch_ours(GEN sel, GEN s, GEN a, GEN lam, long prec)
     859             : {
     860          14 :   GEN selsm0 = gel(sel, 1), selsm1 = gel(sel, 2), selsm2 = gel(sel, 3);
     861          14 :   GEN lin_grid = gel(sel, 4), v0, v1, v2;
     862          14 :   long numpoles = itos(gel(sel, 5));
     863          14 :   GEN worker = snm_closure(is_entry("_serh_worker"),
     864             :                            mkvec4(NULL, NULL, NULL, NULL));
     865          14 :   v0 = add3(series_h0l(worker, m_n0(selsm0), s, a, lam, prec),
     866             :             series_residues_h0l(numpoles, selsm0, s, a, lam, prec),
     867             :             integral_h0l(lin_grid, selsm0, s, a, lam, prec));
     868          14 :   v1 = add3(series_h1(worker, m_n0(selsm1), s, a, lam, prec),
     869             :             series_residues_h1(numpoles, selsm1, s, a, lam, prec),
     870             :             integral_h12(lin_grid, selsm1, s, a, lam, prec));
     871          14 :   v2 = addsub(series_h2(worker, m_n0(selsm2), s, a, lam, prec),
     872             :             series_residues_h2(numpoles, selsm2, s, a, lam, prec),
     873             :             integral_h12(lin_grid, selsm2, s, a, lam, prec));
     874          14 :   return add3(v0, v1, v2);
     875             : }
     876             : 
     877             : static GEN
     878           0 : RZlerch_easy(GEN s, GEN a, GEN lam, long prec)
     879             : {
     880           0 :   pari_sp av = avma;
     881             :   GEN z, y, N;
     882           0 :   long B = prec2nbits(prec), LD = LOWDEFAULTPREC;
     883           0 :   N = gdiv(gmulsg(B + 5, mplog2(LD)), gmul(Pi2n(1, LD), imag_i(lam)));
     884           0 :   if (gexpo(N) > 40) pari_err_IMPL("precision too large in lerchzeta");
     885           0 :   N = gceil(N); prec += EXTRAPREC64;
     886           0 :   z = typ(lam) == t_INT ? gen_1 : gexp(gmul(PiI2(prec), lam), prec);
     887           0 :   y = parsum(gen_0, N, snm_closure(is_entry("_serh_worker"),
     888             :                                    mkvec4(z, a, gneg(s), stoi(prec))));
     889           0 :   return gerepilecopy(av, gprec_wtrunc(y, prec));
     890             : }
     891             : 
     892             : static GEN
     893          14 : mygfrac(GEN z)
     894           0 : { return typ(z) == t_COMPLEX ? mkcomplex(gfrac(real_i(z)), imag_i(z))
     895          14 :                              : gfrac(z); }
     896             : 
     897             : static GEN
     898          42 : lerchlarge(GEN s, GEN a, GEN lam, GEN al, GEN numpoles, long prec)
     899             : {
     900          42 :   pari_sp av = avma;
     901          42 :   GEN val, sel, imlam = imag_i(lam);
     902             :   long prec2;
     903          42 :   switch(gsigne(imlam))
     904             :   {
     905           0 :     case -1: pari_err_IMPL("imag(lam) < 0");
     906           0 :     case  1: if (gexpo(imlam) >= -16) return RZlerch_easy(s, a, lam, prec);
     907             :   }
     908          42 :   if (gcmpgs(real_i(a), 1) < 0)
     909             :   {
     910          14 :     GEN P = gexp(gmul(PiI2(prec), lam), prec);
     911          14 :     GEN L = lerchlarge(s, gaddgs(a, 1), lam, al, numpoles, prec);
     912          14 :     return gerepileupto(av, gadd(gpow(a, gneg(s), prec), gmul(P, L)));
     913             :   }
     914          28 :   if (gcmpgs(real_i(a), 2) >= 0)
     915             :   {
     916           0 :     GEN L, P = gexp(gneg(gmul(PiI2(prec), lam)), prec);
     917           0 :     a = gsubgs(a, 1); L = lerchlarge(s, a, lam, al, numpoles, prec);
     918           0 :     return gerepileupto(av, gmul(P, gsub(L, gpow(a, gneg(s), prec))));
     919             :   }
     920          28 :   if (gsigne(imag_i(s)) > 0)
     921             :   {
     922             :     GEN L;
     923          14 :     lam = mygfrac(gneg(gconj(lam)));
     924          14 :     L = lerchlarge(gconj(s), a, lam, al, numpoles, prec);
     925          14 :     return gerepileupto(av, gconj(L));
     926             :   }
     927          14 :   prec2 = prec + EXTRAPREC64;
     928          14 :   a = gprec_w(a, prec2);
     929          14 :   s = gprec_w(s, prec2);
     930          14 :   lam = gprec_w(lam, prec2);
     931          14 :   sel = RZLERinit(s, a, lam, al, numpoles, prec2);
     932          14 :   val = lerch_ours(sel, s, a, lam, prec2);
     933          14 :   return gerepilecopy(av, gprec_wtrunc(val, prec));
     934             : }
     935             : 
     936             : GEN
     937           7 : zetahurwitzlarge(GEN s, GEN a, long prec)
     938           7 : { return lerchlarge(s, a, gen_0, gen_1, gen_1, prec); }
     939             : 
     940             : GEN
     941           7 : lerchzetalarge(GEN s, GEN a, GEN lam, long prec)
     942           7 : { return lerchlarge(s, a, lam, gen_1, gen_1, prec); }

Generated by: LCOV version 1.16