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.16.2 lcov report (development 29115-f22e516b23) Lines: 297 504 58.9 %
Date: 2024-04-27 08:08:04 Functions: 46 73 63.0 %
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         854 : mycharinit(GEN C, long bit)
      41             : {
      42             :   GEN L, LVC, LE, LGA;
      43         854 :   long F = 0, i, j, lc = lg(C), prec;
      44             : 
      45         854 :   bit += 64; prec = nbits2prec(bit);
      46         854 :   L = cgetg(lc, t_VEC);
      47         854 :   LE = cgetg(lc, t_VECSMALL);
      48         854 :   LGA = cgetg(lc, t_VEC);
      49        1728 :   for (i = 1; i < lc; i++)
      50             :   {
      51         874 :     GEN gv, ga, gm, ro, ldata = gel(C, i);
      52             :     long e;
      53         874 :     if (is_linit(ldata)) ldata = linit_get_ldata(ldata);
      54         874 :     gv = ldata_get_gammavec(ldata); e = itou(gel(gv, 1));
      55         874 :     gm = ldata_get_conductor(ldata);
      56         874 :     ro = ldata_get_rootno(ldata);
      57         874 :     if (isintzero(ro)) ro = lfunrootno(ldata, bit);
      58         874 :     ga = gmul(ro, gsqrt(gm, prec)); if (e) ga = mulcxI(ga);
      59         874 :     gel(LGA, i) = ga;
      60         874 :     LE[i] = e;
      61         874 :     if (i == 1) F = itos(gm); /* constant */
      62         874 :     gel(L, i) = lfunan(ldata, F, prec);
      63             :   }
      64         854 :   if (lc == 2 && is_vec_t(typ(gmael(L,1,1))))
      65             :   { /* multichar */
      66           4 :     LGA = gel(LGA,1); lc = lg(LGA);
      67           4 :     LVC = gel(L,1);
      68           4 :     LE = const_vecsmall(lc-1, LE[1]); /* FIXME: can handle mixed values */
      69             :   }
      70             :   else
      71             :   {
      72         850 :     LVC = cgetg(F + 1, t_VEC);
      73        3949 :     for (j = 1; j <= F; j++)
      74             :     {
      75        3099 :       GEN v = cgetg(lc, t_VEC);
      76        6362 :       for (i = 1; i < lc; i++) gel(v, i) = gmael(L, i, j);
      77        3099 :       gel(LVC, j) = v;
      78             :     }
      79             :   }
      80         854 :   return mkvec4(LVC, stoi(F), LE, LGA);
      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      609259 : mycall(GEN VC, long n)
      87             : {
      88      609259 :   long F = lg(VC) - 1;
      89      609259 :   GEN z = n <= F ? gel(VC, n) : gel(VC, ((n - 1) % F) + 1);
      90      609259 :   return gequal0(z)? NULL: z;
      91             : }
      92             : 
      93      186396 : static GEN get_chivec(GEN VCALL) { return gel(VCALL, 1); }
      94      630774 : static long get_modulus(GEN VCALL) { return itos(gel(VCALL, 2)); }
      95      185361 : static GEN get_signat(GEN VCALL) { return gel(VCALL, 3); }
      96          21 : static GEN get_gauss(GEN VCALL) { return gel(VCALL, 4); }
      97             : 
      98             : /* (-1)^A[i] * conj(B[i]) */
      99             : static GEN
     100           7 : gnegconj(GEN A, GEN B)
     101             : {
     102           7 :   long i, l = lg(A);
     103           7 :   GEN W = cgetg(l, t_VEC);
     104          14 :   for (i = 1; i < l; i++)
     105           7 :   { GEN b = gconj(gel(B,i)); gel(W,i) = A[i]? gneg(b): b; }
     106           7 :   return W;
     107             : }
     108             : /* g(conj(CHI)) */
     109             : static GEN
     110           7 : gaussconj(GEN VCALL)
     111           7 : { return gnegconj(get_signat(VCALL), get_gauss(VCALL)); }
     112             : 
     113             : static GEN
     114           7 : myinitconj(GEN VCALL)
     115             : {
     116           7 :   GEN CONJ = shallowcopy(VCALL);
     117           7 :   gel(CONJ, 1) = gconj(get_chivec(VCALL));
     118           7 :   gel(CONJ, 4) = gaussconj(VCALL); return CONJ;
     119             : }
     120             : 
     121             : /********************************************************************/
     122             : /*                          Driver Program                          */
     123             : /********************************************************************/
     124             : 
     125             : /* assume |Im(s)| >> 1, in particular s is not a negative integer */
     126             : static GEN
     127          14 : applyfuneq(GEN gau, GEN s, GEN z, long odd, long q, long bitprec)
     128             : {
     129             :   GEN t, S;
     130             :   long prec;
     131          14 :   if (!gequal0(s)) bitprec += maxss(gexpo(s), 0);
     132          14 :   prec = nbits2prec(bitprec);
     133          14 :   if (odd) gau = mulcxmI(gau);
     134          14 :   S = gmul(Pi2n(-1, prec), gsubgs(s, odd));
     135          14 :   t = ginv(gmul2n(gmul(gcos(S, prec), ggamma(s, prec)), 1));
     136          14 :   t = gmul(gpow(gdivgs(Pi2n(1, prec), q), s, prec), t);
     137          14 :   return gmul(gmul(gau, t), z);
     138             : }
     139             : 
     140             : static GEN RZchi(GEN VCALL, GEN s, long prec);
     141             : 
     142             : /* VCALL already initialized */
     143             : static GEN
     144         854 : lfunlarge_char(GEN VCALL, GEN s, long bitprec)
     145             : {
     146         854 :   pari_sp av = avma;
     147             :   GEN sig, tau, z;
     148         854 :   long funeq = 0, ts = typ(s), stau, flconj, q;
     149         854 :   if (!is_real_t(ts) && ts != t_COMPLEX) pari_err_TYPE("lfunlarge_char", s);
     150         854 :   sig = real_i(s); tau = imag_i(s);
     151         854 :   if (gexpo(tau) < 1) pari_err_DOMAIN("lfun","im(s)", "<", gen_2, tau);
     152         854 :   stau = gsigne(tau);
     153         854 :   if (stau < 0) { tau = gneg(tau); VCALL = myinitconj(VCALL); }
     154         854 :   if (gcmp(sig, ghalf) < 0) { funeq = 1; sig = gsubsg(1, sig); }
     155         854 :   flconj = ((stau > 0 && funeq) || (stau < 0 && !funeq));
     156         854 :   q = get_modulus(VCALL); bitprec += gexpo(stoi(q));
     157         854 :   z = RZchi(VCALL, mkcomplex(sig, tau), nbits2prec(bitprec));
     158         818 :   if (flconj) z = gconj(z);
     159         818 :   if (funeq)
     160             :   {
     161          14 :     GEN odd = get_signat(VCALL), gau = get_gauss(VCALL), Vz;
     162          14 :     long lC = lg(gau), j;
     163          14 :     Vz = cgetg(lC, t_VEC);
     164          28 :     for (j = 1; j < lC; j++)
     165          14 :       gel(Vz,j) = applyfuneq(gel(gau,j), s, gel(z,j), odd[j], q, bitprec);
     166          14 :     z = Vz;
     167             :   }
     168         818 :   return gerepilecopy(av, z);
     169             : }
     170             : 
     171             : static GEN
     172          14 : lfungetchars(GEN pol)
     173             : {
     174          14 :   GEN L, F, v, bnf = Buchall(pol_x(1), 0, LOWDEFAULTPREC);
     175             :   GEN w, condall, bnr;
     176             :   long i, l, lc;
     177          14 :   condall = rnfconductor(bnf, pol); bnr = gel(condall, 2);
     178          14 :   L = bnrchar(bnr, gel(condall, 3), NULL); lc = lg(L);
     179          14 :   F = cgetg(lc, t_VEC);
     180          77 :   for (i = 1; i < lc; i++)
     181             :   {
     182          63 :     GEN chi = gel(L, i), cond = bnrconductor_raw(bnr, chi);
     183          63 :     gel(F, i) = gcoeff(gel(cond,1), 1, 1);
     184             :   }
     185          14 :   w = vec_equiv(F); l = lg(w); v = cgetg(l, t_COL);
     186          42 :   for (i = 1; i < l; i++)
     187             :   {
     188          28 :     GEN wi = gel(w, i), vi;
     189          28 :     long j, li = lg(wi);
     190          28 :     gel(v,i) = vi = cgetg(li, t_VEC);
     191          28 :     if (li == 2 && equali1(gel(F, wi[1]))) /* common conductor is 1 */
     192          14 :       gel(vi,1) = lfunmisc_to_ldata_shallow(gen_1);
     193             :     else
     194             :     {
     195          63 :       for (j = 1; j < li; j++)
     196          49 :         gel(vi,j) = lfunmisc_to_ldata_shallow(mkvec2(bnr, gel(L, wi[j])));
     197             :     }
     198             :   }
     199          14 :   return v;
     200             : }
     201             : 
     202             : /********************************************************************/
     203             : /*            NEW RS IMPLEMENTATION FROM SANDEEP TYAGI              */
     204             : /********************************************************************/
     205             : /* See arXiv:2203.02509v2 */
     206             : 
     207        4143 : static long m_n0(GEN sel) { return itos(gel(sel, 1)); }
     208      443997 : static GEN m_r0(GEN sel) { return gel(sel, 2); }
     209         854 : static GEN m_al(GEN sel) { return gel(sel, 3); }
     210      443997 : static GEN m_aleps(GEN sel) { return gel(sel, 4); }
     211        3289 : static GEN m_h(GEN sel) { return gel(sel, 5); }
     212         846 : static GEN m_lin(GEN sel) { return gel(sel, 6); }
     213         846 : static GEN m_np(GEN sel) { return gel(sel, 7); }
     214      184299 : static GEN m_pz(GEN sel) { return gel(sel, 8); }
     215             : 
     216             : static GEN
     217        2443 : phi_hat(GEN x, long prec)
     218             : {
     219             :   GEN y;
     220        2443 :   if (signe(imag_i(x)) > 0)
     221         786 :     y = gsubsg(1, gexp(gneg(gmul(PiI2(prec), x)), prec));
     222             :   else
     223        1657 :     y = gsubgs(gexp(gmul(PiI2(prec), x), prec), 1);
     224        2443 :   return ginv(y);
     225             : }
     226             : 
     227             : static GEN
     228        2443 : phi_hat_h0(GEN sel, long k, long prec)
     229             : {
     230        2443 :   GEN t = gdiv(gsubsg(m_n0(sel) + k, m_r0(sel)), m_aleps(sel));
     231        2443 :   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      604452 : mul_addsub(GEN A, GEN a, GEN b, GEN E)
     237             : {
     238      604452 :   long i, l = lg(E);
     239      604452 :   GEN v = cgetg(l, t_VEC);
     240     1257492 :   for (i = 1; i < l; i++)
     241      653040 :     gel(v,i) = gmul(gel(A,i), E[i]? gsub(a, b): gadd(a, b));
     242      604452 :   return v;
     243             : }
     244             : 
     245             : static GEN
     246      184299 : wd(GEN VCALL, GEN pmd, GEN x, GEN PZ, long prec)
     247             : {
     248      184299 :   GEN VC = get_chivec(VCALL), E = get_signat(VCALL), y = NULL;
     249      184299 :   long md = get_modulus(VCALL), k;
     250      788751 :   for (k = 1; k <= (md-1) / 2; k++)
     251             :   {
     252      604452 :     GEN xc = mycall(VC, k);
     253      604452 :     if (xc)
     254             :     {
     255      604452 :       GEN p1 = gmul(xc, gel(PZ, Fl_sqr(k, 2 * md) + 1));
     256      604452 :       GEN p2 = gmul(pmd, gsubgs(x, k)), p3 = gmul(pmd, gaddgs(x, k));
     257      604452 :       p2 = odd(md)? ginv(gsin(p2, prec)): gcotan(p2, prec);
     258      604452 :       p3 = odd(md)? ginv(gsin(p3, prec)): gcotan(p3, prec);
     259      604452 :       p1 = mul_addsub(p1, p2, p3, E);
     260      604452 :       y = y ? gadd(y, p1) : p1;
     261             :     }
     262             :   }
     263      184299 :   return mulcxmI(gdivgs(y, 2*md));
     264             : }
     265             : 
     266             : static GEN
     267         854 : series_h0(long n0, GEN s, GEN VCALL, long fl, long prec)
     268             : {
     269         854 :   GEN C = get_modulus(VCALL) == 1? NULL: get_chivec(VCALL);
     270         854 :   GEN R = pardirpowerssumfun(C, n0, gneg(s), fl, prec);
     271         818 :   if (fl)
     272             :   {
     273          28 :     if (C) return R;
     274          28 :     return mkvec2(mkvec(gel(R,1)), mkvec(gel(R,2)));
     275             :   }
     276         790 :   return C? gel(R,1): mkvec(gel(R,1));
     277             : }
     278             : 
     279             : static GEN
     280         846 : series_residues_h0(GEN sel, GEN s, GEN VCALL, long prec)
     281             : {
     282         846 :   long n0 = m_n0(sel), num_of_poles = itos(m_np(sel)), k;
     283         846 :   GEN val = gen_0, VC = get_chivec(VCALL);
     284        3356 :   for (k = maxss(1 - num_of_poles, 1 - n0); k <= 1 + num_of_poles; k++)
     285             :   {
     286        2510 :     GEN nk = mycall(VC, n0 + k); /* n0 + k > 0 */
     287        2510 :     if (nk) val = gadd(val, gmul(gmul(phi_hat_h0(sel, k, prec), nk),
     288             :                                  gpow(stoi(n0 + k), gneg(s), prec)));
     289             :   }
     290         846 :   return val;
     291             : }
     292             : 
     293             : static GEN
     294      441554 : integrand_h0(GEN sel, GEN s, GEN VCALL, GEN x, long prec)
     295             : {
     296      441554 :   long md = get_modulus(VCALL);
     297      441554 :   GEN r0 = m_r0(sel), aleps = m_aleps(sel), zn, p1;
     298      441554 :   GEN pmd = divru(mppi(prec), md), ix = ginv(x);
     299      441554 :   zn = gadd(r0, gdivgs(gmul(aleps, gsub(x, ix)), 2));
     300      441554 :   p1 = gmul(expIxy(pmd, gsqr(zn), prec),
     301             :             gmul(gpow(zn, gneg(s), prec), gmul(aleps, gadd(x, ix))));
     302      441554 :   if (md == 1)
     303      257255 :     p1 = gdiv(mkvec(mulcxI(p1)), gmul2n(gsin(gmul(pmd, zn), prec), 2));
     304             :   else
     305      184299 :     p1 = gdivgs(gmul(p1, wd(VCALL, pmd, zn, m_pz(sel), prec)), -2);
     306      441554 :   return p1;
     307             : }
     308             : 
     309             : static GEN
     310         846 : integral_h0(GEN sel, GEN s, GEN VCALL, long prec)
     311             : {
     312         846 :   GEN lin_grid = m_lin(sel), S = gen_0;
     313         846 :   pari_sp av = avma;
     314         846 :   long j, l = lg(lin_grid);
     315      411050 :   for (j = 1; j < l; j++)
     316             :   {
     317      410204 :     S = gadd(S, integrand_h0(sel, s, VCALL, gel(lin_grid, j), prec));
     318      410204 :     if ((j & 0xff) == 0) S = gerepileupto(av, S);
     319             :   }
     320         846 :   return gerepileupto(av, gmul(m_h(sel), S));
     321             : }
     322             : 
     323             : /* log |x| */
     324             : static GEN
     325       31350 : mylog(GEN x, long prec)
     326             : {
     327       31350 :   if (gequal0(x)) return gneg(powis(stoi(10), 20)); /* FIXME ! */
     328       31285 :   switch(typ(x))
     329             :   {
     330       31285 :     case t_COMPLEX: return gmul2n(glog(cxnorm(x), prec), -1);
     331           0 :     case t_REAL: break;
     332           0 :     default: x = gtofp(x, prec);
     333             :   }
     334           0 :   return logr_abs(x);
     335             : }
     336             : 
     337             : struct fun_q_t { GEN sel, s, VCALL, B; };
     338             : static GEN
     339       31350 : fun_q(void *E, GEN x)
     340             : {
     341       31350 :   struct fun_q_t *S = (struct fun_q_t *)E;
     342       31350 :   long prec = DEFAULTPREC;
     343       31350 :   GEN z = integrand_h0(S->sel, S->s, S->VCALL, gexp(x, prec), prec);
     344       31350 :   if (typ(z) == t_VEC) z = vecsum(z);
     345       31350 :   return addrr(S->B, mylog(z, prec));
     346             : }
     347             : static GEN
     348        1708 : brent_q(void *E, GEN (*f)(void*,GEN), GEN q_low, GEN q_hi)
     349             : {
     350        1708 :   GEN low_val = f(E, q_low), high_val = f(E, q_hi);
     351        1708 :   if (gsigne(low_val) * gsigne(high_val) >= 0) return NULL;
     352        1643 :   return zbrent(E, f, q_low, q_hi, LOWDEFAULTPREC);
     353             : }
     354             : static GEN
     355         854 : findq(void *E, GEN (*f)(void*,GEN), GEN lq, long B)
     356             : {
     357         854 :   GEN q_low, q_hi, q_right, q_left, q_est = gasinh(lq, LOWDEFAULTPREC);
     358         854 :   q_low = gdivgs(gmulsg(4, q_est), 5);
     359         854 :   q_hi = gdivgs(gmulsg(3, q_est), 2);
     360         854 :   q_right = brent_q(E, f, q_low, q_hi); if (!q_right) q_right = q_est;
     361         854 :   q_left = brent_q(E, f, gneg(q_low), gneg(q_hi)); if (!q_left) q_left = q_est;
     362         854 :   return bitprecision0(gmax(q_right, q_left), B);
     363             : }
     364             : 
     365             : static GEN
     366         854 : set_q_value(GEN sel, GEN s, GEN VCALL, long prec)
     367             : {
     368             :   struct fun_q_t E;
     369         854 :   GEN al = m_al(sel), lq;
     370         854 :   long md = get_modulus(VCALL), LD = DEFAULTPREC;
     371         854 :   E.sel = sel; E.s = s; E.VCALL = VCALL, E.B = mulur(prec, mplog2(LD));
     372         854 :   lq = gdiv(gsqrt(gdiv(gmulsg(md, E.B), Pi2n(1, LD)), LD), al);
     373         854 :   return findq((void*)&E, &fun_q, lq, prec);
     374             : }
     375             : 
     376             : static GEN
     377         854 : setlin_grid_exp(GEN h, long m, long prec)
     378             : {
     379         854 :   GEN w, vex = gpowers(gexp(h, prec), (m - 1)/2);
     380             :   long i;
     381         854 :   w = cgetg(m+1, t_VEC); gel(w, (m + 1)/2) = gen_1;
     382      210685 :   for (i = (m + 3)/2; i <= m; i++)
     383             :   {
     384      209831 :     GEN t1 = gel(vex, i - ((m - 1)/2));
     385      209831 :     gel(w, i) = t1; gel(w, (m + 1) - i) = ginv(t1);
     386             :   }
     387         854 :   return w;
     388             : }
     389             : 
     390             : static long
     391         854 : get_m(GEN q, long prec)
     392             : {
     393         854 :   GEN t = divrr(mulur(4 * prec2nbits(prec), mplog2(prec)), sqrr(mppi(prec)));
     394         854 :   return 2*itos(gfloor(mulrr(q, t))) + 1;
     395             : }
     396             : 
     397             : static GEN
     398         854 : RZinit(GEN s, GEN VCALL, GEN num_of_poles, long prec)
     399             : {
     400             :   GEN sel, al, aleps, n0, r0, q, h, PZ;
     401         854 :   long md = get_modulus(VCALL), m;
     402         854 :   al = gcmpgs(gabs(imag_i(s), prec), 100) < 0 ? ginv(stoi(4)) : gen_1;
     403         854 :   r0 = gsqrt(gdiv(gmulgs(s, md), PiI2(prec)), prec);
     404         854 :   n0 = gfloor(gsub(real_i(r0), imag_i(r0)));
     405         854 :   aleps = gmul(al, gexp(PiI2n(-2, prec), prec));
     406         854 :   PZ = gpowers(gexp(gdivgs(PiI2n(0, prec), -md), prec), 2*md);
     407         854 :   sel = mkvecn(8, n0, r0, al, aleps, NULL, NULL, NULL, PZ);
     408         854 :   q = set_q_value(sel, s, VCALL, prec);
     409         854 :   m = get_m(q, prec);
     410         854 :   gel(sel,5) = h = divru(q, (m - 1) >> 1);
     411         854 :   gel(sel,6) = setlin_grid_exp(h, m, prec);
     412         854 :   gel(sel,7) = num_of_poles;
     413         854 :   return sel;
     414             : }
     415             : 
     416             : /* fl = 0: compute only for s; fl = 1 compute for s and 1-conj(s)
     417             :    and put second result in *ptaux; fl = -1 take *ptaux as serh0. */
     418             : static GEN
     419         882 : total_value(GEN sel, GEN s, GEN VCALL, GEN *ptaux, long fl, long prec)
     420             : {
     421             :   GEN serh0, serh0aux;
     422         882 :   if (fl == -1) serh0 = *ptaux;
     423             :   else
     424             :   {
     425         854 :     serh0aux = series_h0(m_n0(sel), s, VCALL, fl, prec);
     426         818 :     if (fl == 0) serh0 = serh0aux;
     427          28 :     else { serh0 = gel(serh0aux, 1); *ptaux = gel(serh0aux, 2); }
     428             :   }
     429         846 :   return gadd(integral_h0(sel, s, VCALL, prec),
     430             :               gsub(serh0, series_residues_h0(sel, s, VCALL, prec)));
     431             : }
     432             : 
     433             : static GEN
     434         414 : xpquo_one(GEN s, GEN cs, GEN ga, long odd, long md, long prec)
     435             : {
     436         414 :   GEN rho, a = odd? gen_1: gen_0, z = divsr(md, mppi(prec));
     437         414 :   rho = gmul(gdiv(gpow(gen_I(), gdivgs(gneg(a), 2), prec), gsqrt(ga, prec)),
     438             :              gpow(stoi(md), ginv(stoi(4)), prec));
     439         414 :   return gmul(gdiv(gconj(gmul(rho, gpow(z, gdivgs(cs, 2), prec))),
     440             :                    gmul(rho, gpow(z, gdivgs(s, 2), prec))),
     441             :               gexp(gsub(gconj(glngamma(gdivgs(gadd(cs, a), 2), prec)),
     442             :                         glngamma(gdivgs(gadd(s, a), 2), prec)), prec));
     443             : }
     444             : 
     445             : static GEN
     446         854 : xpquo(GEN s, GEN VCALL, long prec)
     447             : {
     448         854 :   long md = get_modulus(VCALL), n, lve, i;
     449         854 :   GEN cd = NULL, z, pz, cs, VC = get_chivec(VCALL), ve, R;
     450         854 :   if (!gequal0(s)) prec = nbits2prec(prec2nbits(prec) + maxss(gexpo(s), 0));
     451         854 :   z = gexp(gdivgs(PiI2(prec), -md), prec);
     452         854 :   if (md == 1)
     453         464 :     return gmul(gpow(mppi(prec), gsub(s, ghalf), prec),
     454             :                 gexp(gsub(glngamma(gdivgs(gsubsg(1, s), 2), prec),
     455             :                           glngamma(gdivgs(s, 2), prec)), prec));
     456         390 :   pz = gpowers(z, md - 1);
     457        2687 :   for (n = 1; n < md; n++)
     458             :   {
     459        2297 :     GEN xn = mycall(VC, n);
     460        2297 :     if (xn)
     461             :     {
     462        2262 :       GEN tmp = gmul(xn, gel(pz, n + 1));
     463        2262 :       cd = cd ? gadd(cd, tmp) : tmp;
     464             :     }
     465             :   }
     466         390 :   cs = gsubsg(1, gconj(s));
     467         390 :   ve = get_signat(VCALL); lve = lg(ve); R = cgetg(lve, t_VEC);
     468         804 :   for (i = 1; i < lve; i++)
     469         414 :     gel(R, i) = xpquo_one(s, cs, gel(cd, i), ve[i], md, prec);
     470         390 :   if (lve == 2) R = gel(R, 1);
     471         390 :   return R;
     472             : }
     473             : 
     474             : static GEN
     475         854 : dirichlet_ours(GEN s, GEN VCALL, long prec)
     476             : {
     477         854 :   GEN sel = RZinit(s, VCALL, gen_1, prec);
     478         854 :   GEN cs, sum1, sum2, aux, pre_fac = xpquo(s, VCALL, prec);
     479         854 :   if (gequal(real_i(s), ghalf))
     480         826 :   { sum1 = total_value(sel, s, VCALL, NULL, 0, prec); sum2 = sum1; }
     481             :   else
     482             :   {
     483          28 :     sum1 = total_value(sel, s, VCALL, &aux, 1, prec);
     484          28 :     cs = gsubsg(1, gconj(s));
     485          28 :     sum2 = total_value(sel, cs, VCALL, &aux, -1, prec);
     486             :   }
     487         818 :   return gadd(sum1, vecmul(pre_fac, gconj(sum2)));
     488             : }
     489             : 
     490             : /* assume |Im(s)| > 2^-bitprec */
     491             : static GEN
     492         854 : RZchi(GEN VCALL, GEN s, long prec)
     493             : {
     494         854 :   long prec2 = prec + EXTRAPRECWORD;
     495         854 :   return gprec_wtrunc(dirichlet_ours(gprec_w(s, prec2), VCALL, prec2), prec);
     496             : }
     497             : 
     498             : /********************************************************************/
     499             : /*                         Utility Functions                        */
     500             : /********************************************************************/
     501             : /* lam = 0, return L(s); else Lambda(s) */
     502             : static GEN
     503         854 : lfuncharall(GEN VCALL, GEN s, long lam, long bitprec)
     504             : {
     505         854 :   GEN ve, P, Q, R, z = lfunlarge_char(VCALL, s, bitprec);
     506             :   long l, i, q, prec;
     507         818 :   if (!lam) return z;
     508         651 :   ve = get_signat(VCALL); l = lg(ve);
     509         651 :   q = get_modulus(VCALL); prec = nbits2prec(bitprec);
     510         651 :   R = cgetg(l, t_VEC);
     511         651 :   Q = divur(q, mppi(prec));
     512         651 :   P = (q == 1 || zv_equal0(ve))? NULL: gsqrt(utoipos(q), prec);
     513        1302 :   for (i = 1; i < l; i++)
     514             :   {
     515         651 :     GEN se = gmul2n(gaddgs(s, ve[i]), -1), r;
     516         651 :     if (lam == 1)
     517             :     {
     518           0 :       r = gmul(gpow(Q, se, prec), ggamma(se, prec));
     519           0 :       if (P && ve[i]) r = gdiv(r, P);
     520             :     }
     521             :     else
     522             :     {
     523         651 :       r = gadd(gmul(se, glog(Q, prec)), glngamma(se, prec));
     524         651 :       if (P && ve[i]) r = gsub(r, glog(P, prec));
     525             :     }
     526         651 :     gel(R, i) = r;
     527             :   }
     528         651 :   return lam == 1 ? vecmul(R, z) : gadd(R, glog(z, prec));
     529             : }
     530             : 
     531             : static GEN
     532          28 : lfunlargeall_from_chars(GEN v, GEN s, long lam, long bit)
     533             : {
     534          28 :   long i, l = lg(v);
     535          72 :   for (i = 1; i < l; i++)
     536             :   {
     537          50 :     GEN w = mycharinit(gel(v, i), bit), L = lfuncharall(w, s, lam, bit);
     538          44 :     gel(v, i) = lam==-1 ? vecsum(L): vecprod(L);
     539             :   }
     540          22 :   return lam==-1 ? vecsum(v): vecprod(v);
     541             : }
     542             : static GEN
     543         832 : lfunlargeall(GEN ldata, GEN s, long lam, long bit)
     544             : {
     545             :   GEN w, an;
     546         832 :   if (lg(ldata) == 2)
     547             :   { /* HACK: ldata[1] a t_DESC_PRODUCT from lfunabelianrelinit / Q */
     548          14 :     GEN v = lfunprod_get_fact(linit_get_tech(gel(ldata,1)));
     549             :     long i, l;
     550          14 :     v = shallowcopy(gel(v,1)); l = lg(v);
     551          42 :     for (i = 1; i < l; i++) gel(v,i) = mkvec(gel(v,i));
     552          14 :     return lfunlargeall_from_chars(v, s, lam, bit);
     553             :   }
     554         818 :   an = gel(ldata_get_an(ldata), 2);
     555         818 :   switch(ldata_get_type(ldata))
     556             :   {
     557          14 :     case t_LFUN_NF:
     558             :     {
     559          14 :       GEN v = lfungetchars(nf_get_pol(an));
     560          14 :       return lfunlargeall_from_chars(v, s, lam, bit);
     561             :     }
     562           4 :     case t_LFUN_CHIGEN:
     563             :     {
     564           4 :       GEN chi = gmael(an, 2, 2);
     565           4 :       if (lg(chi) > 1 && is_vec_t(typ(gel(chi,1))))
     566             :       { /* multi char */
     567           4 :         w = mycharinit(mkcol(ldata), bit);
     568           4 :         return lfuncharall(w, s, lam, bit);
     569             :       }
     570             :     }
     571             :     default: /* single char */
     572         800 :       w = mycharinit(mkcol(ldata), bit);
     573         800 :       return gel(lfuncharall(w, s, lam, bit), 1);
     574             :   }
     575             : }
     576             : 
     577             : GEN
     578         183 : lfunlarge(GEN CHI, GEN s, long bit)
     579         183 : { return lfunlargeall(CHI, s, 0, bit); }
     580             : 
     581             : GEN
     582           0 : lfunlambdalarge(GEN CHI, GEN s, long bit)
     583           0 : { return lfunlargeall(CHI, s, 1, bit); }
     584             : 
     585             : GEN
     586         649 : lfunloglambdalarge(GEN CHI, GEN s, long bit)
     587         649 : { return lfunlargeall(CHI, s, -1, bit); }
     588             : 
     589             : /********************************************************************/
     590             : /*           LERCH RS IMPLEMENTATION FROM SANDEEP TYAGI             */
     591             : /********************************************************************/
     592             : 
     593             : static GEN
     594           0 : double_exp_residue_pos_h(GEN selsm, long k, long ind, long prec)
     595             : {
     596           0 :   long nk = itos(gel(selsm, 1)) + k;
     597           0 :   GEN r = gel(selsm, 2), ale = gel(selsm, 3), aor = gel(selsm, 4);
     598           0 :   GEN h = gel(selsm, 5), t = gen_0;
     599           0 :   switch(ind)
     600             :   {
     601           0 :     case 0: t = gaddsg(nk, aor); break;
     602           0 :     case 1: t = gneg(gaddsg(nk, aor)); break;
     603           0 :     case 2: t = gsubsg(nk, aor); break;
     604             :   }
     605           0 :   return gdiv(gasinh(gdiv(gsub(t, r), ale), prec), h);
     606             : }
     607             : 
     608             : static GEN
     609           0 : phi_hat_h(GEN selsm, long m, long ind, long prec)
     610           0 : { return phi_hat(double_exp_residue_pos_h(selsm, m, ind, prec), prec); }
     611             : 
     612             : static long
     613           0 : myex(GEN is)
     614           0 : { return gequal0(is) ? 0 : maxss(0, 2 + gexpo(is)); }
     615             : 
     616             : static GEN
     617           0 : gaminus(GEN s, long prec)
     618             : {
     619           0 :   GEN is = imag_i(s), tmp;
     620           0 :   long B = prec2nbits(prec);
     621           0 :   if (gcmpgs(is, -5*B) < 0) return gen_0;
     622           0 :   prec = nbits2prec(B + myex(is));
     623           0 :   tmp = gexp(gsub(glngamma(s, prec), gmul(PiI2n(-1, prec), s)), prec);
     624           0 :   return bitprecision0(tmp, B);
     625             : }
     626             : 
     627             : static GEN
     628           0 : gaplus(GEN s, long prec)
     629             : {
     630           0 :   GEN is = imag_i(s), tmp;
     631           0 :   long B = prec2nbits(prec);
     632           0 :   if (gcmpgs(is, 5*B) > 0) return gen_0;
     633           0 :   prec = nbits2prec(B + myex(is));
     634           0 :   tmp = gexp(gadd(glngamma(s, prec), gmul(PiI2n(-1, prec), s)), prec);
     635           0 :   return bitprecision0(tmp, B);
     636             : }
     637             : 
     638             : GEN
     639           0 : serh_worker(GEN gk, GEN z, GEN a, GEN ns, GEN gprec)
     640             : {
     641           0 :   long k = itos(gk);
     642           0 :   return gmul(gpowgs(z, k), gpow(gaddsg(k, a), ns, itos(gprec)));
     643             : }
     644             : 
     645             : static GEN
     646           0 : allparsums(long ini, long fin, GEN z, GEN a, GEN ns, long prec)
     647             : {
     648           0 :   return parsum(stoi(ini), stoi(fin), strtoclosure("_serh_worker", 4, z, a, ns, stoi(prec)));
     649             : }
     650             : 
     651             : static GEN
     652           0 : series_h0l(long n0, GEN s, GEN a, GEN lam, long prec)
     653             : {
     654           0 :   GEN z = typ(lam) == t_INT ? gen_1 : gexp(gmul(PiI2(prec), lam), prec);
     655           0 :   return allparsums(0, n0, z, a, gneg(s), prec);
     656             : }
     657             : 
     658             : static GEN
     659           0 : series_h1(long n1, GEN s, GEN a, GEN lam, long prec)
     660             : {
     661           0 :   GEN sum1, pre_factor, z, sn = gsubgs(s, 1);
     662           0 :   GEN ini = gequal0(lam) ? gen_1 : gen_0;
     663           0 :   pre_factor = gaplus(gneg(sn), prec);
     664           0 :   if (gequal0(pre_factor)) return gen_0;
     665           0 :   pre_factor = gmul(gmul(pre_factor, gexp(gneg(gmul(PiI2(prec), gmul(a, lam))), prec)), gpow(Pi2n(1, prec), sn, prec));
     666           0 :   z = typ(a) == t_INT ? gen_1 : gexp(gneg(gmul(PiI2(prec), a)), prec);
     667           0 :   sum1 = allparsums(itos(ini), n1 - 1, z, lam, sn, prec);
     668           0 :   return gmul(pre_factor, sum1);
     669             : }
     670             : 
     671             : static GEN
     672           0 : series_h2(long n2, GEN s, GEN a, GEN lam, long prec)
     673             : {
     674           0 :   GEN sum2, pre_factor, z, sn = gsubgs(s, 1);
     675           0 :   pre_factor = gaminus(gneg(sn), prec);
     676           0 :   if (gequal0(pre_factor)) return gen_0;
     677           0 :   pre_factor = gmul(gmul(pre_factor, gexp(gneg(gmul(PiI2(prec), gmul(a, lam))), prec)), gpow(Pi2n(1, prec), sn, prec));
     678           0 :   z = typ(a) == t_INT ? gen_1 : gexp(gmul(PiI2(prec), a), prec);
     679           0 :   sum2 = allparsums(1, n2, z, gneg(lam), sn, prec);
     680           0 :   return gmul(pre_factor, sum2);
     681             : }
     682             : 
     683             : static GEN
     684           0 : series_residues_h0l(long num_of_poles, GEN selsm0, GEN s, GEN a, GEN lam, long prec)
     685             : {
     686           0 :   GEN val = gen_0, ra = real_i(a);
     687           0 :   long n0 = m_n0(selsm0), k;
     688           0 :   for (k = -num_of_poles + 1; k <= num_of_poles; k++)
     689             :   {
     690           0 :     if (gsigne(gaddsg(n0 + k, ra)) > 0)
     691           0 :       val = gadd(val, gmul(gmul(phi_hat_h(selsm0, k, 0, prec), gexp(gmul(PiI2(prec), gmulgs(lam, n0 + k)), prec)), gpow(gaddsg(n0 + k, a), gneg(s), prec)));
     692             :   }
     693           0 :   return val;
     694             : }
     695             : 
     696             : static GEN
     697           0 : series_residues_h1(long num_of_poles, GEN selsm1, GEN s, GEN a, GEN lam, long prec)
     698             : {
     699           0 :   GEN val = gen_0, rlam = real_i(lam), pre_factor, sn = gsubgs(s, 1);
     700           0 :   long n1 = m_n0(selsm1), k;
     701           0 :   pre_factor = gaplus(gneg(sn), prec);
     702           0 :   if (gequal0(pre_factor)) return gen_0;
     703           0 :   pre_factor = gmul(gmul(pre_factor, gexp(gneg(gmul(PiI2(prec), gmul(a, lam))), prec)), gpow(Pi2n(1, prec), sn, prec));
     704           0 :   for (k = -num_of_poles; k <= num_of_poles - 1; k++)
     705             :   {
     706           0 :     if (gsigne(gaddsg(n1 + k, rlam)) > 0)
     707           0 :       val = gadd(val, gmul(gmul(phi_hat_h(selsm1, k, 1, prec), gexp(gneg(gmul(PiI2(prec), gmulgs(a, n1 + k))), prec)), gpow(gaddsg(n1 + k, lam), sn, prec)));
     708             :   }
     709           0 :   return gmul(pre_factor, val);
     710             : }
     711             : 
     712             : static GEN
     713           0 : series_residues_h2(long num_of_poles, GEN selsm2, GEN s, GEN a, GEN lam, long prec)
     714             : {
     715           0 :   GEN val = gen_0, rlam = real_i(lam), pre_factor, sn = gsubgs(s, 1);
     716           0 :   long n2 = m_n0(selsm2), k;
     717           0 :   pre_factor = gaminus(gneg(sn), prec);
     718           0 :   if (gequal0(pre_factor)) return gen_0;
     719           0 :   pre_factor = gmul(gmul(pre_factor, gexp(gneg(gmul(PiI2(prec), gmul(a, lam))), prec)), gpow(Pi2n(1, prec), sn, prec));
     720           0 :   for (k = -num_of_poles + 1; k <= num_of_poles; k++)
     721             :   {
     722           0 :     if (gsigne(gsubsg(n2 + k, rlam)) > 0)
     723           0 :       val = gsub(val, gmul(gmul(phi_hat_h(selsm2, k, 2, prec), gexp(gmul(PiI2(prec), gmulgs(a, n2 + k)), prec)), gpow(gsubsg(n2 + k, lam), sn, prec)));
     724             :   }
     725           0 :   return gmul(pre_factor, val);
     726             : }
     727             : 
     728             : static GEN
     729           0 : integrand_h0l(GEN selsm0, GEN s, GEN alam, GEN x, long prec)
     730             : {
     731           0 :   GEN r0 = gel(selsm0, 2), ale = gel(selsm0, 3), a = gel(selsm0, 4);
     732           0 :   GEN ix = ginv(x), zn = gadd(r0, gmul2n(gmul(ale, gsub(x, ix)), -1));
     733           0 :   GEN pii = PiI2n(0, prec), den, num;
     734           0 :   den = gsub(gexp(gmul(pii, gsub(zn, gmul2n(a, 1))), prec), gexp(gneg(gmul(pii, zn)), prec));
     735           0 :   num = gexp(gmul(gmul(pii, zn), gsub(gmul2n(alam, 1), zn)), prec);
     736           0 :   num = gmul(gmul(gmul(num, ale), gmul2n(gadd(x, ix), -1)), gpow(zn, gneg(s), prec));
     737           0 :   return gdiv(num, den);
     738             : }
     739             : 
     740             : static GEN
     741           0 : integrand_h12(GEN selsm1, GEN s, GEN alam, GEN x, long prec)
     742             : {
     743           0 :   GEN r1 = gel(selsm1, 2), ale = gel(selsm1, 3), lam = gel(selsm1, 4);
     744           0 :   GEN ix = ginv(x), zn = gadd(r1, gmul2n(gmul(ale, gsub(x, ix)), -1));
     745           0 :   GEN pii = PiI2n(0, prec), den, num, y;
     746           0 :   den = gsub(gexp(gmul(pii, gadd(zn, gmul2n(lam, 1))), prec), gexp(gneg(gmul(pii, zn)), prec));
     747           0 :   num = gexp(gmul(gmul(pii, zn), gadd(gmul2n(alam, 1), zn)), prec);
     748           0 :   num = gmul(gmul(gmul(num, ale), gmul2n(gadd(x, ix), -1)), gpow(zn, gsubgs(s, 1), prec));
     749           0 :   y = gdiv(num, den);
     750           0 :   if (gcmp(garg(zn, prec), Pi2n(-2, prec)) > 0)
     751           0 :     y = gmul(y, gexp(gmul(PiI2(prec), gsubsg(1, s)), prec));
     752           0 :   return y;
     753             : }
     754             : 
     755             : static GEN
     756           0 : integral_h0l(GEN lin_grid, GEN selsm0, GEN s, GEN a, GEN lam, long prec)
     757             : {
     758           0 :   GEN alam = gadd(a, lam), S = gen_0;
     759           0 :   pari_sp av = avma;
     760           0 :   long j, l = lg(lin_grid);
     761           0 :   for (j = 1; j < l; j++)
     762             :   {
     763           0 :     S = gadd(S, integrand_h0l(selsm0, s, alam, gel(lin_grid, j), prec));
     764           0 :     if ((j & 0xff) == 0) S = gerepileupto(av, S);
     765             :   }
     766           0 :   S = gmul(m_h(selsm0), S);
     767           0 :   return gmul(S, gexp(gneg(gmul(PiI2n(0, prec), gmul(a, gaddsg(1, gadd(a, gmul2n(lam, 1)))))), prec));
     768             : }
     769             : 
     770             : /* do not forget a minus sign for index 2 */
     771             : static GEN
     772           0 : integral_h12(GEN lin_grid, GEN selsm1, GEN s, GEN a, GEN lam, long prec)
     773             : {
     774           0 :   GEN alam = gadd(a, lam), S = gen_0, ga = gaminus(gsubsg(1, s), prec);
     775           0 :   pari_sp av = avma;
     776           0 :   long j, l = lg(lin_grid);
     777           0 :   if (gequal0(ga)) return S;
     778           0 :   for (j = 1; j < l; j++)
     779             :   {
     780           0 :     S = gadd(S, integrand_h12(selsm1, s, alam, gel(lin_grid, j), prec));
     781           0 :     if ((j & 0xff) == 0) S = gerepileupto(av, S);
     782             :   }
     783           0 :   if (gequal0(S)) return gen_0;
     784           0 :   S = gmul(m_h(selsm1), S);
     785           0 :   return gmul(gmul(gmul(S, ga), gexp(gmul(PiI2n(0, prec), gmul(lam, gaddgs(lam, 1))), prec)), gpow(Pi2n(1, prec), gsubgs(s, 1), prec));
     786             : }
     787             : 
     788             : struct _fun_q0_t { GEN sel, s, alam, B; };
     789             : static GEN
     790           0 : _fun_q0(void *E, GEN x)
     791             : {
     792           0 :   struct _fun_q0_t *S = (struct _fun_q0_t*)E;
     793           0 :   GEN z = integrand_h0l(S->sel, S->s, S->alam, x, DEFAULTPREC);
     794           0 :   return addrr(S->B, mylog(z, DEFAULTPREC));
     795             : }
     796             : static GEN
     797           0 : _fun_q12(void *E, GEN x)
     798             : {
     799           0 :   struct _fun_q0_t *S = (struct _fun_q0_t*)E;
     800           0 :   GEN z = integrand_h12(S->sel, S->s, S->alam, x, DEFAULTPREC);
     801           0 :   return addrr(S->B, mylog(z, DEFAULTPREC));
     802             : }
     803             : 
     804             : static GEN
     805           0 : RZLERinit(GEN s, GEN a, GEN lam, GEN al, GEN num_of_poles, long prec)
     806             : {
     807             :   GEN eps, r0, r1, r2, h, lin_grid, q, q0, q1, q2, sel0, sel1, sel2, lq;
     808           0 :   GEN pinv = ginv(PiI2(prec)), c = gmul2n(gadd(a, lam), -1), n0, n1, n2, c2;
     809           0 :   long m, LD = DEFAULTPREC;
     810             :   struct _fun_q0_t E;
     811           0 :   if (!al || gequal0(al))
     812           0 :     al = gcmpgs(gabs(imag_i(s), prec), 100) < 0 ? ginv(stoi(4)) : gen_1;
     813           0 :   c2 = gsub(gsqr(c), gmul(s, pinv));
     814           0 :   r0 = gadd(c, gsqrt(c2, prec));
     815           0 :   r1 = gsqrt(gadd(c2, pinv), prec);
     816           0 :   r2 = gsub(r1, c);
     817           0 :   r1 = gneg(gadd(r1, c));
     818           0 :   n0 = gfloor(gsub(gadd(real_i(r0), imag_i(r0)), a));
     819           0 :   n1 = gneg(gfloor(gadd(gsub(real_i(r1), imag_i(r1)), lam)));
     820           0 :   n2 = gfloor(gadd(gsub(real_i(r2), imag_i(r2)), lam));
     821             : 
     822           0 :   eps = gexp(PiI2n(-2, prec), prec);
     823           0 :   E.s = s; E.alam = gadd(a, lam), E.B = mulur(prec, mplog2(prec));
     824           0 :   lq = gmul(gsqrt(gmul(gdivsg(prec, Pi2n(1, LD)), mplog2(LD)), LD), al);
     825           0 :   E.sel = sel0 = mkvec5(n0, r0, gdiv(al, eps), a, gen_0);
     826           0 :   q0 = findq(&E, &_fun_q0, lq, prec);
     827             : 
     828           0 :   if (!gequal1(al)) lq = gdiv(lq, gsqr(al));
     829           0 :   E.sel = sel1 = mkvec5(n1, r1, gmul(al, eps), lam, gen_0);
     830           0 :   q1 = findq(&E, &_fun_q12, lq, prec);
     831           0 :   E.sel = sel2 = mkvec5(n2, r2, gmul(al, eps), lam, gen_0);
     832           0 :   q2 = findq(&E, &_fun_q12, lq, prec);
     833           0 :   q = vecmax(mkvec3(q0, q1, q2));
     834           0 :   m = get_m(q, prec);
     835           0 :   gel(sel0, 5) = gel(sel1, 5) = gel(sel2, 5) = h = divru(q, (m-1) >> 1);
     836           0 :   lin_grid = setlin_grid_exp(h, m, prec);
     837           0 :   if (!num_of_poles) num_of_poles = gen_1;
     838           0 :   return mkvec5(sel0, sel1, sel2, lin_grid, num_of_poles);
     839             : }
     840             : 
     841             : static GEN
     842           0 : lerch_ours(GEN sel, GEN s, GEN a, GEN lam, long prec)
     843             : {
     844           0 :   GEN selsm0 = gel(sel, 1), selsm1 = gel(sel, 2), selsm2 = gel(sel, 3);
     845           0 :   GEN lin_grid = gel(sel, 4), val_h0, val_h1, val_h2;
     846             :   GEN serandres_h0, serandres_h1, serandres_h2;
     847           0 :   long num_of_poles = itos(gel(sel, 5));
     848           0 :   serandres_h0 = gadd(series_h0l(m_n0(selsm0), s, a, lam, prec), series_residues_h0l(num_of_poles, selsm0, s, a, lam, prec));
     849           0 :   val_h0 = gadd(integral_h0l(lin_grid, selsm0, s, a, lam, prec), serandres_h0);
     850           0 :   serandres_h1 = gadd(series_h1(m_n0(selsm1), s, a, lam, prec), series_residues_h1(num_of_poles, selsm1, s, a, lam, prec));
     851           0 :   val_h1 = gadd(integral_h12(lin_grid, selsm1, s, a, lam, prec), serandres_h1);  serandres_h2 = gadd(series_h2(m_n0(selsm2), s, a, lam, prec), series_residues_h2(num_of_poles, selsm2, s, a, lam, prec));
     852           0 :   val_h2 = gadd(gneg(integral_h12(lin_grid, selsm2, s, a, lam, prec)), serandres_h2);
     853           0 :   return gadd(gadd(val_h0, val_h1), val_h2);
     854             : }
     855             : 
     856             : GEN
     857           0 : serhlong_worker(GEN gk, GEN z, GEN a, GEN ns, GEN gprec)
     858             : {
     859           0 :   long k = itos(gk);
     860           0 :   return gmul(gpowgs(z, k), gpow(gaddsg(k, a), ns, itos(gprec)));
     861             : }
     862             : 
     863             : static GEN
     864           0 : RZlerch_easy(GEN s, GEN a, GEN lam, long prec)
     865             : {
     866           0 :   pari_sp ltop = avma;
     867             :   GEN z, y, gnlim;
     868           0 :   long B = prec2nbits(prec), nlim, LD = LOWDEFAULTPREC;
     869           0 :   gnlim = gceil(gdiv(gmulsg(B + 5, mplog2(LD)), gmul(Pi2n(1, prec), imag_i(lam))));
     870           0 :   if (gexpo(gnlim) > 40)
     871           0 :     pari_err(e_MISC, "imag(lam) too small for RZlerch_easy");
     872           0 :   nlim = itos(gnlim); prec += EXTRAPRECWORD;
     873           0 :   z = typ(lam) == t_INT ? gen_1 : gexp(gmul(PiI2(prec), lam), prec);
     874           0 :   if (nlim < 10000000) y = allparsums(0, nlim, z, a, gneg(s), prec);
     875             :   else
     876           0 :     y = parsum(gen_0, gnlim, strtoclosure("_serhlong_worker", 4, z, a, gneg(s), stoi(prec)));
     877           0 :   return gerepileupto(ltop, bitprecision0(y, B));
     878             : }
     879             : 
     880             : static GEN
     881           0 : lerchlarge(GEN s, GEN a, GEN lam, GEN al, GEN num_of_poles, long prec)
     882             : {
     883           0 :   pari_sp ltop = avma;
     884             :   GEN val, sel;
     885           0 :   long B = prec2nbits(prec), NB = B + EXTRAPRECWORD, sl = gsigne(imag_i(lam));
     886           0 :   if (sl < 0) pari_err_IMPL("imag(lam) < 0");
     887           0 :   if (sl > 0) return RZlerch_easy(s, a, lam, prec);
     888           0 :   if (gcmpgs(real_i(a), 1) < 0)
     889           0 :     return gerepileupto(ltop, gadd(gpow(a, gneg(s), prec), gmul(gexp(gmul(PiI2(prec), lam), prec), lerchlarge(s, gaddgs(a, 1), lam, al, num_of_poles, prec))));
     890           0 :   if (gcmpgs(real_i(a), 2) >= 0)
     891           0 :     return gerepileupto(ltop, gmul(gexp(gneg(gmul(PiI2(prec), lam)), prec), gsub(lerchlarge(s, gsubgs(a, 1), lam, al, num_of_poles, prec), gpow(gsubgs(a, 1), gneg(s), prec))));
     892           0 :   if (gsigne(imag_i(s)) > 0)
     893           0 :     return gerepileupto(ltop, gconj(lerchlarge(gconj(s), a, gfrac(gneg(lam)), al, num_of_poles, prec)));
     894           0 :   a = bitprecision0(a, NB);
     895           0 :   s = bitprecision0(s, NB);
     896           0 :   lam = bitprecision0(lam, NB); prec = nbits2prec(NB);
     897           0 :   sel = RZLERinit(s, a, lam, al, num_of_poles, prec);
     898           0 :   val = lerch_ours(sel, s, a, lam, prec);
     899           0 :   return gerepilecopy(ltop, bitprecision0(val, B));
     900             : }
     901             : 
     902             : GEN
     903           0 : zetahurwitzlarge(GEN s, GEN a, long prec)
     904           0 : { return lerchlarge(s, a, gen_0, gen_1, gen_1, prec); }
     905             : 
     906             : GEN
     907           0 : lerchzetalarge(GEN s, GEN a, GEN lam, long prec)
     908           0 : { return lerchlarge(s, a, lam, gen_1, gen_1, prec); }

Generated by: LCOV version 1.14