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 - lfun.c (source / functions) Hit Total Coverage
Test: PARI/GP v2.12.1 lcov report (development 24215-a79de5b25) Lines: 1392 1447 96.2 %
Date: 2019-08-24 05:50:50 Functions: 142 143 99.3 %
Legend: Lines: hit not hit

          Line data    Source code
       1             : /* Copyright (C) 2015  The PARI group.
       2             : 
       3             : This file is part of the PARI/GP package.
       4             : 
       5             : PARI/GP is free software; you can redistribute it and/or modify it under the
       6             : terms of the GNU General Public License as published by the Free Software
       7             : Foundation. It is distributed in the hope that it will be useful, but WITHOUT
       8             : ANY WARRANTY WHATSOEVER.
       9             : 
      10             : Check the License for details. You should have received a copy of it, along
      11             : with the package; see the file 'COPYING'. If not, write to the Free Software
      12             : Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. */
      13             : 
      14             : /********************************************************************/
      15             : /**                                                                **/
      16             : /**                       L-functions                              **/
      17             : /**                                                                **/
      18             : /********************************************************************/
      19             : 
      20             : #include "pari.h"
      21             : #include "paripriv.h"
      22             : 
      23             : /*******************************************************************/
      24             : /*  Accessors                                                      */
      25             : /*******************************************************************/
      26             : 
      27             : static GEN
      28       11889 : mysercoeff(GEN x, long n)
      29             : {
      30       11889 :   long N = n - valp(x);
      31       11889 :   return (N < 0)? gen_0: gel(x, N+2);
      32             : }
      33             : 
      34             : long
      35        5642 : ldata_get_type(GEN ldata) { return mael3(ldata, 1, 1, 1); }
      36             : 
      37             : GEN
      38       14616 : ldata_get_an(GEN ldata) { return gel(ldata, 1); }
      39             : 
      40             : GEN
      41       32879 : ldata_get_dual(GEN ldata) { return gel(ldata, 2); }
      42             : 
      43             : long
      44        1713 : ldata_isreal(GEN ldata) { return isintzero(gel(ldata, 2)); }
      45             : 
      46             : GEN
      47      189584 : ldata_get_gammavec(GEN ldata) { return gel(ldata, 3); }
      48             : 
      49             : long
      50       11779 : ldata_get_degree(GEN ldata) { return lg(gel(ldata, 3))-1; }
      51             : 
      52             : GEN
      53       88738 : ldata_get_k(GEN ldata)
      54             : {
      55       88738 :   GEN w = gel(ldata,4);
      56       88738 :   if (typ(w) == t_VEC) w = gel(w,1);
      57       88738 :   return w;
      58             : }
      59             : /* a_n = O(n^{k1 + epsilon}) */
      60             : static double
      61       63391 : ldata_get_k1(GEN ldata)
      62             : {
      63       63391 :   GEN w = gel(ldata,4);
      64             :   double k;
      65       63391 :   if (typ(w) == t_VEC) return gtodouble(gel(w,2));
      66             :   /* by default, assume that k1 = k-1 and even (k-1)/2 for entire functions */
      67       63048 :   k = gtodouble(w);
      68       63048 :   return ldata_get_residue(ldata)? k-1: (k-1)/2.;
      69             : }
      70             : 
      71             : GEN
      72      133667 : ldata_get_conductor(GEN ldata) { return gel(ldata, 5); }
      73             : 
      74             : GEN
      75       45629 : ldata_get_rootno(GEN ldata) { return gel(ldata, 6); }
      76             : 
      77             : GEN
      78      114102 : ldata_get_residue(GEN ldata) { return lg(ldata) == 7 ? NULL: gel(ldata, 7); }
      79             : 
      80             : long
      81       82201 : linit_get_type(GEN linit) { return mael(linit, 1, 1); }
      82             : 
      83             : GEN
      84      120725 : linit_get_ldata(GEN linit) { return gel(linit, 2); }
      85             : 
      86             : GEN
      87      152965 : linit_get_tech(GEN linit) { return gel(linit, 3); }
      88             : 
      89             : long
      90      121408 : is_linit(GEN data)
      91             : {
      92      220643 :   return lg(data) == 4 && typ(data) == t_VEC
      93      220640 :                        && typ(gel(data, 1)) == t_VECSMALL;
      94             : }
      95             : 
      96             : GEN
      97       18903 : lfun_get_step(GEN tech) { return gmael(tech, 2, 1);}
      98             : 
      99             : GEN
     100       18903 : lfun_get_pol(GEN tech) { return gmael(tech, 2, 2);}
     101             : 
     102             : GEN
     103        4247 : lfun_get_Residue(GEN tech) { return gmael(tech, 2, 3);}
     104             : 
     105             : GEN
     106       29975 : lfun_get_k2(GEN tech) { return gmael(tech, 3, 1);}
     107             : 
     108             : GEN
     109       11513 : lfun_get_w2(GEN tech) { return gmael(tech, 3, 2);}
     110             : 
     111             : GEN
     112       11513 : lfun_get_expot(GEN tech) { return gmael(tech, 3, 3);}
     113             : 
     114             : GEN
     115        4039 : lfun_get_factgammavec(GEN tech) { return gmael(tech, 3, 4); }
     116             : 
     117             : /* Handle complex Vga whose sum is real */
     118             : static GEN
     119       68578 : sumVga(GEN Vga) { return real_i(vecsum(Vga)); }
     120             : /* sum_i max (Im v[i],0) */
     121             : static double
     122       15552 : sumVgaimpos(GEN v)
     123             : {
     124       15552 :   double d = 0.;
     125       15552 :   long i, l = lg(v);
     126       40822 :   for (i = 1; i < l; i++)
     127             :   {
     128       25270 :     GEN c = imag_i(gel(v,i));
     129       25270 :     if (gsigne(c) > 0) d += gtodouble(c);
     130             :   }
     131       15552 :   return d;
     132             : }
     133             : 
     134             : static long
     135       14148 : vgaell(GEN Vga)
     136             : {
     137             :   GEN c;
     138       14148 :   long d = lg(Vga)-1;
     139       14148 :   if (d != 2) return 0;
     140        9185 :   c = gsub(gel(Vga,1), gel(Vga,2));
     141        9185 :   return gequal1(c) || gequalm1(c);
     142             : }
     143             : static long
     144       47867 : vgaeasytheta(GEN Vga) { return lg(Vga)-1 == 1 || vgaell(Vga); }
     145             : /* return b(n) := a(n) * n^c, when vgaeasytheta(Vga) is set */
     146             : static GEN
     147       10087 : antwist(GEN an, GEN Vga, long prec)
     148             : {
     149             :   long l, i;
     150       10087 :   GEN b, c = vecmin(Vga);
     151       10087 :   if (gequal0(c)) return an;
     152        1064 :   l = lg(an); b = cgetg(l, t_VEC);
     153        1064 :   if (gequal1(c))
     154             :   {
     155         693 :     if (typ(an) == t_VECSMALL)
     156         245 :       for (i = 1; i < l; i++) gel(b,i) = mulss(an[i], i);
     157             :     else
     158         448 :       for (i = 1; i < l; i++) gel(b,i) = gmulgs(gel(an,i), i);
     159             :   }
     160             :   else
     161             :   {
     162         371 :     GEN v = vecpowug(l-1, c, prec);
     163         371 :     if (typ(an) == t_VECSMALL)
     164           0 :       for (i = 1; i < l; i++) gel(b,i) = gmulsg(an[i], gel(v,i));
     165             :     else
     166         371 :       for (i = 1; i < l; i++) gel(b,i) = gmul(gel(an,i), gel(v,i));
     167             :   }
     168        1064 :   return b;
     169             : }
     170             : 
     171             : static GEN
     172        6391 : theta_dual(GEN theta, GEN bn)
     173             : {
     174        6391 :   if (typ(bn)==t_INT) return NULL;
     175             :   else
     176             :   {
     177         119 :     GEN thetad = shallowcopy(theta), ldata = linit_get_ldata(theta);
     178         119 :     GEN Vga = ldata_get_gammavec(ldata);
     179         119 :     GEN tech = shallowcopy(linit_get_tech(theta));
     180         119 :     GEN an = theta_get_an(tech);
     181         119 :     long prec = nbits2prec(theta_get_bitprec(tech));
     182         119 :     GEN vb = ldata_vecan(bn, lg(an)-1, prec);
     183         119 :     if (!theta_get_m(tech) && vgaeasytheta(Vga)) vb = antwist(vb, Vga, prec);
     184         119 :     gel(tech,1) = vb;
     185         119 :     gel(thetad,3) = tech; return thetad;
     186             :   }
     187             : }
     188             : 
     189             : static GEN
     190       33041 : domain_get_dom(GEN domain)  { return gel(domain,1); }
     191             : static long
     192       14516 : domain_get_der(GEN domain)  { return mael2(domain, 2, 1); }
     193             : static long
     194       19227 : domain_get_bitprec(GEN domain)  { return mael2(domain, 2, 2); }
     195             : GEN
     196       33468 : lfun_get_domain(GEN tech) { return gel(tech,1); }
     197             : long
     198          49 : lfun_get_bitprec(GEN tech){ return domain_get_bitprec(lfun_get_domain(tech)); }
     199             : GEN
     200           0 : lfun_get_dom(GEN tech) { return domain_get_dom(lfun_get_domain(tech)); }
     201             : 
     202             : GEN
     203        1797 : lfunprod_get_fact(GEN tech)  { return gel(tech, 2); }
     204             : 
     205             : GEN
     206       40993 : theta_get_an(GEN tdata)      { return gel(tdata, 1);}
     207             : GEN
     208        5425 : theta_get_K(GEN tdata)       { return gel(tdata, 2);}
     209             : GEN
     210        4755 : theta_get_R(GEN tdata)       { return gel(tdata, 3);}
     211             : long
     212       54127 : theta_get_bitprec(GEN tdata) { return itos(gel(tdata, 4));}
     213             : long
     214       81335 : theta_get_m(GEN tdata)       { return itos(gel(tdata, 5));}
     215             : GEN
     216       43261 : theta_get_tdom(GEN tdata)    { return gel(tdata, 6);}
     217             : GEN
     218       45354 : theta_get_isqrtN(GEN tdata)  { return gel(tdata, 7);}
     219             : 
     220             : /*******************************************************************/
     221             : /*  Helper functions related to Gamma products                     */
     222             : /*******************************************************************/
     223             : 
     224             : /* return -itos(s) >= 0 if s is (approximately) equal to a non-positive
     225             :  * integer, and -1 otherwise */
     226             : static long
     227       15519 : isnegint(GEN s)
     228             : {
     229       15519 :   GEN r = ground(real_i(s));
     230       15519 :   if (signe(r) <= 0 && gequal(s, r)) return -itos(r);
     231       15477 :   return -1;
     232             : }
     233             : 
     234             : /* pi^(-s/2) Gamma(s/2) */
     235             : static GEN
     236       10997 : gamma_R(GEN s, long prec)
     237             : {
     238       10997 :   GEN s2 = gdivgs(s, 2), pi = mppi(prec);
     239       10997 :   long ms = isnegint(s2);
     240       10997 :   if (ms >= 0)
     241             :   {
     242          42 :     GEN pr = gmul(powru(pi, ms), gdivsg(odd(ms)? -2: 2, mpfact(ms)));
     243          42 :     GEN S = scalarser(pr, 0, 1);
     244          42 :     setvalp(S,-1); return S;
     245             :   }
     246       10955 :   return gdiv(ggamma(s2,prec), gpow(pi,s2,prec));
     247             : }
     248             : 
     249             : /* gamma_R(s)gamma_R(s+1) = 2 (2pi)^(-s) Gamma(s) */
     250             : static GEN
     251        4522 : gamma_C(GEN s, long prec)
     252             : {
     253        4522 :   GEN pi2 = Pi2n(1,prec);
     254        4522 :   long ms = isnegint(s);
     255        4522 :   if (ms >= 0)
     256             :   {
     257           0 :     GEN pr = gmul(powrs(pi2, ms), gdivsg(odd(ms)? -2: 2, mpfact(ms)));
     258           0 :     GEN S = scalarser(pr, 0, 1);
     259           0 :     setvalp(S,-1); return S;
     260             :   }
     261        4522 :   return gmul2n(gdiv(ggamma(s,prec), gpow(pi2,s,prec)), 1);
     262             : }
     263             : 
     264             : static GEN
     265        1099 : gammafrac(GEN r, long d)
     266             : {
     267        1099 :   GEN pr, a = gmul2n(r, -1);
     268        1099 :   GEN polj = cgetg(labs(d)+1, t_COL);
     269        1099 :   long i, v=0;
     270        1099 :   if (d > 0)
     271           0 :     for (i = 1; i <= d; ++i)
     272           0 :       gel(polj, i) = deg1pol_shallow(ghalf, gaddgs(a, i-1), v);
     273             :   else
     274        2198 :     for (i = 1; i <= -d; ++i)
     275        1099 :       gel(polj, i) = deg1pol_shallow(ghalf, gsubgs(a, i), v);
     276        1099 :   pr = RgV_prod(polj);
     277        1099 :   return d < 0 ? ginv(pr): pr;
     278             : }
     279             : 
     280             : static GEN
     281       15729 : gammafactor(GEN Vga)
     282             : {
     283       15729 :   pari_sp av = avma;
     284       15729 :   long i, m, d = lg(Vga)-1, dr, dc;
     285       15729 :   GEN pol = pol_1(0), pi = gen_0, R = cgetg(d+1,t_VEC);
     286             :   GEN P, F, FR, FC, E, ER, EC;
     287       40264 :   for (i = 1; i <= d; ++i)
     288             :   {
     289       24535 :     GEN a = gel(Vga,i), qr = gdiventres(real_i(a), gen_2);
     290       24535 :     long q = itos(gel(qr,1));
     291       24535 :     gel(R, i) = gadd(gel(qr,2), imag_i(a));
     292       24535 :     if (q)
     293             :     {
     294        1099 :       pol = gmul(pol, gammafrac(gel(R,i), q));
     295        1099 :       pi  = addis(pi, q);
     296             :     }
     297             :   }
     298       15729 :   gen_sort_inplace(R, (void*)cmp_universal, cmp_nodata, &P);
     299       15729 :   F = cgetg(d+1, t_VEC); E = cgetg(d+1, t_VECSMALL);
     300       52577 :   for (i = 1, m = 0; i <= d;)
     301             :   {
     302             :     long k;
     303       21119 :     GEN u = gel(R, i);
     304       24535 :     for(k = i + 1; k <= d; ++k)
     305        8806 :       if (cmp_universal(gel(R, k), u)) break;
     306       21119 :     m++;
     307       21119 :     E[m] = k - i;
     308       21119 :     gel(F, m) = u;
     309       21119 :     i = k;
     310             :   }
     311       15729 :   setlg(F, m+1); setlg(E, m+1);
     312       15729 :   R = cgetg(m+1, t_VEC);
     313       36848 :   for (i = 1; i <= m; i++)
     314             :   {
     315       21119 :     GEN qr = gdiventres(gel(F,i), gen_1);
     316       21119 :     gel(R, i) = mkvec2(gel(qr,2), stoi(E[i]));
     317             :   }
     318       15729 :   gen_sort_inplace(R, (void*)cmp_universal, cmp_nodata, &P);
     319       15729 :   FR = cgetg(m+1, t_VEC); ER = cgetg(m+1, t_VECSMALL);
     320       15729 :   FC = cgetg(m+1, t_VEC); EC = cgetg(m+1, t_VECSMALL);
     321       47936 :   for (i = 1, dr = 1, dc = 1; i <= m;)
     322             :   {
     323       16478 :     if (i==m || cmp_universal(gel(R,i), gel(R,i+1)))
     324             :     {
     325       11837 :       gel(FR, dr) = gel(F, P[i]);
     326       11837 :       ER[dr] = E[P[i]];
     327       11837 :       dr++; i++;
     328             :     } else
     329             :     {
     330        4641 :       if (gequal(gaddgs(gmael(R,i,1), 1), gmael(R,i+1,1)))
     331           0 :         gel(FC, dc) = gel(F, P[i+1]);
     332             :       else
     333        4641 :         gel(FC, dc) = gel(F, P[i]);
     334        4641 :       EC[dc] = E[P[i]];
     335        4641 :       dc++; i+=2;
     336             :     }
     337             :   }
     338       15729 :   setlg(FR, dr); setlg(ER, dr);
     339       15729 :   setlg(FC, dc); setlg(EC, dc);
     340       15729 :   return gerepilecopy(av, mkvec4(pol, pi, mkvec2(FR,ER), mkvec2(FC,EC)));
     341             : }
     342             : 
     343             : static GEN
     344       13755 : deg1ser_shallow(GEN a1, GEN a0, long v, long e)
     345             : {
     346       13755 :   return RgX_to_ser(deg1pol_shallow(a1, a0, v), e+2);
     347             : }
     348             : /*
     349             : To test:
     350             : GR(s)=Pi^-(s/2)*gamma(s/2);
     351             : GC(s)=2*(2*Pi)^-s*gamma(s)
     352             : gam_direct(F,s)=prod(i=1,#F,GR(s+F[i]))
     353             : gam_fact(F,s)=my([P,p,R,C]=gammafactor(F));subst(P,x,s)*Pi^-p*prod(i=1,#R[1],GR(s+R[1][i])^R[2][i])*prod(i=1,#C[1],GC(s+C[1][i])^C[2][i])
     354             : */
     355             : 
     356             : static GEN
     357       15918 : polgammaeval(GEN F, GEN s)
     358             : {
     359       15918 :   GEN r = poleval(F, s);
     360       15918 :   if (typ(s)!=t_SER && gequal0(r))
     361             :   {
     362           0 :     long e = gvaluation(F, deg1pol(gen_1, gneg(s), varn(F)));
     363           0 :     r = poleval(F, deg1ser_shallow(gen_1, s, 0, e+1));
     364             :   }
     365       15918 :   return r;
     366             : }
     367             : 
     368             : static GEN
     369       14728 : fracgammaeval(GEN F, GEN s)
     370             : {
     371       14728 :   if (typ(F)==t_POL)
     372       13538 :     return polgammaeval(F, s);
     373        1190 :   else if (typ(F)==t_RFRAC)
     374        1190 :     return gdiv(polgammaeval(gel(F,1), s), polgammaeval(gel(F,2), s));
     375           0 :   return F;
     376             : }
     377             : 
     378             : static GEN
     379       14728 : gammafactproduct(GEN F, GEN s, long prec)
     380             : {
     381       14728 :   pari_sp av = avma;
     382       14728 :   GEN P = fracgammaeval(gel(F,1), s);
     383       14728 :   GEN p = gpow(mppi(prec),gneg(gel(F,2)), prec), z = gmul(P, p);
     384       14728 :   GEN R = gel(F,3), Rw = gel(R,1), Re=gel(R,2);
     385       14728 :   GEN C = gel(F,4), Cw = gel(C,1), Ce=gel(C,2);
     386       14728 :   long i, lR = lg(Rw), lC = lg(Cw);
     387       25725 :   for (i=1; i< lR; i++)
     388       10997 :     z = gmul(z, gpowgs(gamma_R(gadd(s,gel(Rw, i)), prec), Re[i]));
     389       19250 :   for (i=1; i< lC; i++)
     390        4522 :     z = gmul(z, gpowgs(gamma_C(gadd(s,gel(Cw, i)), prec), Ce[i]));
     391       14728 :   return gerepileupto(av, z);
     392             : }
     393             : 
     394             : static int
     395        4060 : gammaordinary(GEN Vga, GEN s)
     396             : {
     397        4060 :   long i, d = lg(Vga)-1;
     398       10850 :   for (i = 1; i <= d; i++)
     399             :   {
     400        6874 :     GEN z = gadd(s, gel(Vga,i));
     401             :     long e;
     402        6874 :     if (gsigne(z) <= 0) { (void)grndtoi(z, &e); if (e < -10) return 0; }
     403             :   }
     404        3976 :   return 1;
     405             : }
     406             : 
     407             : /* Exponent A of t in asymptotic expansion; K(t) ~ C t^A exp(-pi d t^(2/d)).
     408             :  * suma = vecsum(Vga)*/
     409             : static double
     410       63384 : gammavec_expo(long d, double suma) { return (1 - d + suma) / d; }
     411             : 
     412             : /*******************************************************************/
     413             : /*       First part: computations only involving Theta(t)          */
     414             : /*******************************************************************/
     415             : 
     416             : static void
     417      102685 : get_cone(GEN t, double *r, double *a)
     418             : {
     419      102685 :   const long prec = LOWDEFAULTPREC;
     420      102685 :   if (typ(t) == t_COMPLEX)
     421             :   {
     422       15218 :     t  = gprec_w(t, prec);
     423       15218 :     *r = gtodouble(gabs(t, prec));
     424       15218 :     *a = fabs(gtodouble(garg(t, prec)));
     425             :   }
     426             :   else
     427             :   {
     428       87467 :     *r = fabs(gtodouble(t));
     429       87467 :     *a = 0.;
     430             :   }
     431      102685 :   if (!*r && !*a) pari_err_DOMAIN("lfunthetainit","t","=",gen_0,t);
     432      102678 : }
     433             : /* slightly larger cone than necessary, to avoid round-off problems */
     434             : static void
     435       59424 : get_cone_fuzz(GEN t, double *r, double *a)
     436       59424 : { get_cone(t, r, a); *r -= 1e-10; if (*a) *a += 1e-10; }
     437             : 
     438             : /* Initialization m-th Theta derivative. tdom is either
     439             :  * - [rho,alpha]: assume |t| >= rho and |arg(t)| <= alpha
     440             :  * - a positive real scalar: assume t real, t >= tdom;
     441             :  * - a complex number t: compute at t;
     442             :  * N is the conductor (either the true one from ldata or a guess from
     443             :  * lfunconductor) */
     444             : long
     445       47839 : lfunthetacost(GEN ldata, GEN tdom, long m, long bitprec)
     446             : {
     447       47839 :   pari_sp av = avma;
     448       47839 :   GEN Vga = ldata_get_gammavec(ldata);
     449       47839 :   long d = lg(Vga)-1;
     450       47839 :   long k1 = ldata_get_k1(ldata);
     451       47839 :   double c = d/2., a, A, B, logC, al, rho, T;
     452       47839 :   double N = gtodouble(ldata_get_conductor(ldata));
     453             : 
     454       47839 :   if (!N) pari_err_TYPE("lfunthetaneed [missing conductor]", ldata);
     455       47839 :   if (typ(tdom) == t_VEC && lg(tdom) == 3)
     456             :   {
     457           7 :     rho= gtodouble(gel(tdom,1));
     458           7 :     al = gtodouble(gel(tdom,2));
     459             :   }
     460             :   else
     461       47832 :     get_cone_fuzz(tdom, &rho, &al);
     462       47832 :   A = gammavec_expo(d, gtodouble(sumVga(Vga))); set_avma(av);
     463       47832 :   a = (A+k1+1) + (m-1)/c;
     464       47832 :   if (fabs(a) < 1e-10) a = 0.;
     465       47832 :   logC = c*M_LN2 - log(c)/2;
     466             :   /* +1: fudge factor */
     467       47832 :   B = M_LN2*bitprec+logC+m*log(2*M_PI) + 1 + (k1+1)*log(N)/2 - (k1+m+1)*log(rho);
     468       47832 :   if (al)
     469             :   { /* t = rho e^(i*al), T^(1/c) = Re(t^(1/c)) > 0, T = rho cos^c(al/c) */
     470        7609 :     double z = cos(al/c);
     471        7609 :     T = (d == 2 && typ(tdom) != t_VEC)? gtodouble(real_i(tdom)): rho*pow(z,c);
     472        7609 :     if (z <= 0)
     473           0 :       pari_err_DOMAIN("lfunthetaneed", "arg t", ">", dbltor(c*M_PI/2), tdom);
     474        7609 :     B -= log(z) * (c * (k1+A+1) + m);
     475             :   }
     476             :   else
     477       40223 :     T = rho;
     478       47832 :   return B <= 0? 0: floor(0.9 + dblcoro526(a,c,B) / T * sqrt(N));
     479             : }
     480             : long
     481          14 : lfunthetacost0(GEN L, GEN tdom, long m, long bitprec)
     482             : {
     483             :   long n;
     484          14 :   if (is_linit(L) && linit_get_type(L)==t_LDESC_THETA)
     485           7 :   {
     486           7 :     GEN tech = linit_get_tech(L);
     487           7 :     n = lg(theta_get_an(tech))-1;
     488             :   }
     489             :   else
     490             :   {
     491           7 :     pari_sp av = avma;
     492           7 :     GEN ldata = lfunmisc_to_ldata_shallow(L);
     493           7 :     n = lfunthetacost(ldata, tdom? tdom: gen_1, m, bitprec);
     494           7 :     set_avma(av);
     495             :   }
     496          14 :   return n;
     497             : }
     498             : 
     499             : static long
     500        8022 : fracgammadegree(GEN FVga)
     501        8022 : { GEN F = gel(FVga,1); return (typ(F)==t_RFRAC)? degpol(gel(F,2)): 0; }
     502             : 
     503             : /* Poles of a L-function can be represented in the following ways:
     504             :  * 1) Nothing (ldata has only 6 components, ldata_get_residue = NULL).
     505             :  * 2) a complex number (single pole at s = k with given residue, unknown if 0).
     506             :  * 3) A vector (possibly empty) of 2-component vectors [a, ra], where a is the
     507             :  * pole, ra a t_SER: its Taylor expansion at a. A t_VEC encodes the polar
     508             :  * part of L, a t_COL, the polar part of Lambda */
     509             : 
     510             : /* 'a' a complex number (pole), 'r' the polar part of L at 'a';
     511             :  * return 'R' the polar part of Lambda at 'a' */
     512             : static GEN
     513        6657 : rtoR(GEN a, GEN r, GEN FVga, GEN N, long prec)
     514             : {
     515        6657 :   long v = lg(r)-2;
     516        6657 :   GEN as = deg1ser_shallow(gen_1, a, varn(r), v);
     517        6657 :   GEN Na = gpow(N, gdivgs(as, 2), prec);
     518        6657 :   long d = fracgammadegree(FVga);
     519        6657 :   if (d) as = sertoser(as, v+d); /* make up for a possible loss of accuracy */
     520        6657 :   return gmul(gmul(r, Na), gammafactproduct(FVga, as, prec));
     521             : }
     522             : 
     523             : /* assume r in normalized form: t_VEC of pairs [be,re] */
     524             : GEN
     525        6615 : lfunrtopoles(GEN r)
     526             : {
     527        6615 :   long j, l = lg(r);
     528        6615 :   GEN v = cgetg(l, t_VEC);
     529       13384 :   for (j = 1; j < l; j++)
     530             :   {
     531        6769 :     GEN rj = gel(r,j), a = gel(rj,1);
     532        6769 :     gel(v,j) = a;
     533             :   }
     534        6615 :   gen_sort_inplace(v, (void*)&cmp_universal, cmp_nodata, NULL);
     535        6615 :   return v;
     536             : }
     537             : 
     538             : /* r / x + O(1) */
     539             : static GEN
     540        5439 : simple_pole(GEN r)
     541             : {
     542             :   GEN S;
     543        5439 :   if (isintzero(r)) return gen_0;
     544        5432 :   S = deg1ser_shallow(gen_0, r, 0, 1);
     545        5432 :   setvalp(S, -1); return S;
     546             : }
     547             : static GEN
     548        5845 : normalize_simple_pole(GEN r, GEN k)
     549             : {
     550        5845 :   long tx = typ(r);
     551        5845 :   if (is_vec_t(tx)) return r;
     552        5439 :   if (!is_scalar_t(tx)) pari_err_TYPE("lfunrootres [poles]", r);
     553        5439 :   return mkvec(mkvec2(k, simple_pole(r)));
     554             : }
     555             : /* normalize the description of a polar part */
     556             : static GEN
     557        7504 : normalizepoles(GEN r, GEN k)
     558             : {
     559             :   long iv, j, l;
     560             :   GEN v;
     561        7504 :   if (!is_vec_t(typ(r))) return normalize_simple_pole(r, k);
     562        3150 :   v = cgetg_copy(r, &l);
     563        7455 :   for (j = iv = 1; j < l; j++)
     564             :   {
     565        4305 :     GEN rj = gel(r,j), a = gel(rj,1), ra = gel(rj,2);
     566        4305 :     if (!is_scalar_t(typ(a)) || typ(ra) != t_SER)
     567           0 :       pari_err_TYPE("lfunrootres [poles]",r);
     568        4305 :     if (valp(ra) >= 0) continue;
     569        4305 :     gel(v,iv++) = rj;
     570             :   }
     571        3150 :   setlg(v, iv); return v;
     572             : }
     573             : static int
     574        9163 : residues_known(GEN r)
     575             : {
     576        9163 :   long i, l = lg(r);
     577        9163 :   if (isintzero(r)) return 0;
     578        9037 :   if (!is_vec_t(typ(r))) return 1;
     579       10584 :   for (i = 1; i < l; i++)
     580             :   {
     581        6055 :     GEN ri = gel(r,i);
     582        6055 :     if (!is_vec_t(typ(ri)) || lg(ri)!=3)
     583           0 :       pari_err_TYPE("lfunrootres [poles]",r);
     584        6055 :     if (isintzero(gel(ri, 2))) return 0;
     585             :   }
     586        4529 :   return 1;
     587             : }
     588             : 
     589             : /* Compute R's from r's (r = Taylor devts of L(s), R of Lambda(s)).
     590             :  * 'r/eno' passed to override the one from ldata  */
     591             : static GEN
     592       16912 : lfunrtoR_i(GEN ldata, GEN r, GEN eno, long prec)
     593             : {
     594       16912 :   GEN Vga = ldata_get_gammavec(ldata), N = ldata_get_conductor(ldata);
     595             :   GEN R, vr, FVga;
     596       16912 :   pari_sp av = avma;
     597             :   long lr, j, jR;
     598       16912 :   GEN k = ldata_get_k(ldata);
     599             : 
     600       16912 :   if (!r || isintzero(eno) || !residues_known(r))
     601        9408 :     return gen_0;
     602        7504 :   r = normalizepoles(r, k);
     603        7504 :   if (typ(r) == t_COL) return gerepilecopy(av, r);
     604        6503 :   if (typ(ldata_get_dual(ldata)) != t_INT)
     605           0 :     pari_err(e_MISC,"please give the Taylor development of Lambda");
     606        6503 :   vr = lfunrtopoles(r); lr = lg(vr);
     607        6503 :   FVga = gammafactor(Vga);
     608        6503 :   R = cgetg(2*lr, t_VEC);
     609       13160 :   for (j = jR = 1; j < lr; j++)
     610             :   {
     611        6657 :     GEN rj = gel(r,j), a = gel(rj,1), ra = gel(rj,2);
     612        6657 :     GEN Ra = rtoR(a, ra, FVga, N, prec);
     613        6657 :     GEN b = gsub(k, conj_i(a));
     614        6657 :     if (lg(Ra)-2 < -valp(Ra))
     615           0 :       pari_err(e_MISC,
     616             :         "please give more terms in L function's Taylor development at %Ps", a);
     617        6657 :     gel(R,jR++) = mkvec2(a, Ra);
     618        6657 :     if (!tablesearch(vr, b, (int (*)(GEN,GEN))&cmp_universal))
     619             :     {
     620        6580 :       GEN mX = gneg(pol_x(varn(Ra)));
     621        6580 :       GEN Rb = gmul(eno, gsubst(conj_i(Ra), varn(Ra), mX));
     622        6580 :       gel(R,jR++) = mkvec2(b, Rb);
     623             :     }
     624             :   }
     625        6503 :   setlg(R, jR); return gerepilecopy(av, R);
     626             : }
     627             : static GEN
     628       16828 : lfunrtoR_eno(GEN ldata, GEN eno, long prec)
     629       16828 : { return lfunrtoR_i(ldata, ldata_get_residue(ldata), eno, prec); }
     630             : static GEN
     631       14070 : lfunrtoR(GEN ldata, long prec)
     632       14070 : { return lfunrtoR_eno(ldata, ldata_get_rootno(ldata), prec); }
     633             : 
     634             : /* thetainit using {an: n <= L}; if (m = 0 && easytheta), an2 is an * n^al */
     635             : static GEN
     636       11599 : lfunthetainit0(GEN ldata, GEN tdom, GEN an2, long m,
     637             :     long bitprec, long extrabit)
     638             : {
     639       11599 :   long prec = nbits2prec(bitprec);
     640       11599 :   GEN tech, N = ldata_get_conductor(ldata);
     641       11599 :   GEN Vga = ldata_get_gammavec(ldata);
     642       11599 :   GEN K = gammamellininvinit(Vga, m, bitprec + extrabit);
     643       11599 :   GEN R = lfunrtoR(ldata, prec);
     644       11599 :   if (!tdom) tdom = gen_1;
     645       11599 :   if (typ(tdom) != t_VEC)
     646             :   {
     647             :     double r, a;
     648       11592 :     get_cone_fuzz(tdom, &r, &a);
     649       11592 :     tdom = mkvec2(dbltor(r), a? dbltor(a): gen_0);
     650             :   }
     651       11599 :   tech = mkvecn(7, an2,K,R, stoi(bitprec), stoi(m), tdom,
     652             :                    gsqrt(ginv(N), prec + EXTRAPRECWORD));
     653       11599 :   return mkvec3(mkvecsmall(t_LDESC_THETA), ldata, tech);
     654             : }
     655             : 
     656             : /* tdom: 1) positive real number r, t real, t >= r; or
     657             :  *       2) [r,a], describing the cone |t| >= r, |arg(t)| <= a */
     658             : static GEN
     659        7098 : lfunthetainit_i(GEN data, GEN tdom, long m, long bitprec)
     660             : {
     661        7098 :   GEN ldata = lfunmisc_to_ldata_shallow(data);
     662        7098 :   long L = lfunthetacost(ldata, tdom, m, bitprec), prec = nbits2prec(bitprec);
     663        7091 :   GEN an = ldata_vecan(ldata_get_an(ldata), L, prec);
     664        7091 :   GEN Vga = ldata_get_gammavec(ldata);
     665        7091 :   if (m == 0 && vgaeasytheta(Vga)) an = antwist(an, Vga, prec);
     666        7091 :   return lfunthetainit0(ldata, tdom, an, m, bitprec, 32);
     667             : }
     668             : 
     669             : GEN
     670         259 : lfunthetainit(GEN ldata, GEN tdom, long m, long bitprec)
     671             : {
     672         259 :   pari_sp av = avma;
     673         259 :   GEN S = lfunthetainit_i(ldata, tdom? tdom: gen_1, m, bitprec);
     674         259 :   return gerepilecopy(av, S);
     675             : }
     676             : 
     677             : GEN
     678         840 : lfunan(GEN ldata, long L, long prec)
     679             : {
     680         840 :   pari_sp av = avma;
     681             :   GEN an ;
     682         840 :   ldata = lfunmisc_to_ldata_shallow(ldata);
     683         840 :   an = gerepilecopy(av, ldata_vecan(ldata_get_an(ldata), L, prec));
     684         833 :   if (typ(an) != t_VEC) an = vecsmall_to_vec_inplace(an);
     685         833 :   return an;
     686             : }
     687             : 
     688             : /* [1^B,...,N^B] */
     689             : GEN
     690         231 : vecpowuu(long N, ulong B)
     691             : {
     692             :   GEN v;
     693             :   long p, i;
     694             :   forprime_t T;
     695             : 
     696         231 :   if (B <= 2)
     697             :   {
     698          63 :     if (!B) return const_vec(N,gen_1);
     699          56 :     v = cgetg(N+1, t_VEC); if (N == 0) return v;
     700          56 :     gel(v,1) = gen_1;
     701          56 :     if (B == 1)
     702          42 :       for (i = 2; i <= N; i++) gel(v,i) = utoipos(i);
     703             :     else
     704          14 :       for (i = 2; i <= N; i++) gel(v,i) = sqru(i);
     705          56 :     return v;
     706             :   }
     707         168 :   v = const_vec(N, NULL);
     708         168 :   u_forprime_init(&T, 3, N);
     709        3752 :   while ((p = u_forprime_next(&T)))
     710             :   {
     711             :     long m, pk, oldpk;
     712        3416 :     gel(v,p) = powuu(p, B);
     713        7623 :     for (pk = p, oldpk = p; pk; oldpk = pk, pk = umuluu_le(pk,p,N))
     714             :     {
     715        4207 :       if (pk != p) gel(v,pk) = mulii(gel(v,oldpk), gel(v,p));
     716       18844 :       for (m = N/pk; m > 1; m--)
     717       14637 :         if (gel(v,m) && m%p) gel(v, m*pk) = mulii(gel(v,m), gel(v,pk));
     718             :     }
     719             :   }
     720         168 :   gel(v,1) = gen_1;
     721        6909 :   for (i = 2; i <= N; i+=2)
     722             :   {
     723        6741 :     long vi = vals(i);
     724        6741 :     gel(v,i) = shifti(gel(v,i >> vi), B * vi);
     725             :   }
     726         168 :   return v;
     727             : }
     728             : /* [1^B,...,N^B] */
     729             : GEN
     730        8544 : vecpowug(long N, GEN B, long prec)
     731             : {
     732        8544 :   GEN v = const_vec(N, NULL);
     733        8544 :   long p, eB = gexpo(B);
     734        8544 :   long prec0 = eB < 5? prec: prec + nbits2extraprec(eB);
     735             :   forprime_t T;
     736        8544 :   u_forprime_init(&T, 2, N);
     737        8544 :   gel(v,1) = gen_1;
     738      345198 :   while ((p = u_forprime_next(&T)))
     739             :   {
     740             :     long m, pk, oldpk;
     741      328110 :     gel(v,p) = gpow(utor(p,prec0), B, prec);
     742      328110 :     if (prec0 != prec) gel(v,p) = gprec_wtrunc(gel(v,p), prec);
     743      727432 :     for (pk = p, oldpk = p; pk; oldpk = pk, pk = umuluu_le(pk,p,N))
     744             :     {
     745      399322 :       if (pk != p) gel(v,pk) = gmul(gel(v,oldpk), gel(v,p));
     746     6034872 :       for (m = N/pk; m > 1; m--)
     747     5635550 :         if (gel(v,m) && m%p) gel(v, m*pk) = gmul(gel(v,m), gel(v,pk));
     748             :     }
     749             :   }
     750        8544 :   return v;
     751             : }
     752             : 
     753             : GEN
     754          49 : dirpowers(long n, GEN x, long prec)
     755             : {
     756          49 :   pari_sp av = avma;
     757             :   GEN v;
     758          49 :   if (n <= 0) return cgetg(1, t_VEC);
     759          35 :   if (typ(x) == t_INT && lgefint(x) <= 3 && signe(x) >= 0)
     760           7 :   {
     761          28 :     ulong B = itou(x);
     762          28 :     v = vecpowuu(n, B);
     763          28 :     if (B <= 2) return v;
     764             :   }
     765           7 :   else v = vecpowug(n, x, prec);
     766          14 :   return gerepilecopy(av, v);
     767             : }
     768             : 
     769             : /* return [1^(2/d), 2^(2/d),...,lim^(2/d)] */
     770             : static GEN
     771        5425 : mkvroots(long d, long lim, long prec)
     772             : {
     773        5425 :   if (d <= 4)
     774             :   {
     775        5271 :     GEN v = cgetg(lim+1,t_VEC);
     776             :     long n;
     777        5271 :     switch(d)
     778             :     {
     779             :       case 1:
     780        1883 :         for (n=1; n <= lim; n++) gel(v,n) = sqru(n);
     781        1883 :         return v;
     782             :       case 2:
     783         987 :         for (n=1; n <= lim; n++) gel(v,n) = utoipos(n);
     784         987 :         return v;
     785             :       case 4:
     786        1382 :         for (n=1; n <= lim; n++) gel(v,n) = sqrtr(utor(n, prec));
     787        1382 :         return v;
     788             :     }
     789             :   }
     790        1173 :   return vecpowug(lim, gdivgs(gen_2,d), prec);
     791             : }
     792             : 
     793             : GEN
     794       47209 : lfunthetacheckinit(GEN data, GEN t, long m, long bitprec)
     795             : {
     796       47209 :   if (is_linit(data) && linit_get_type(data)==t_LDESC_THETA)
     797             :   {
     798       43261 :     GEN tdom, thetainit = linit_get_tech(data);
     799       43261 :     long bitprecnew = theta_get_bitprec(thetainit);
     800       43261 :     long m0 = theta_get_m(thetainit);
     801             :     double r, al, rt, alt;
     802       43261 :     if (m0 != m)
     803           0 :       pari_err_DOMAIN("lfuntheta","derivative order","!=", stoi(m),stoi(m0));
     804       43261 :     if (bitprec > bitprecnew) goto INIT;
     805       43261 :     get_cone(t, &rt, &alt);
     806       43261 :     tdom = theta_get_tdom(thetainit);
     807       43261 :     r = rtodbl(gel(tdom,1));
     808       43261 :     al= rtodbl(gel(tdom,2)); if (rt >= r && alt <= al) return data;
     809             :   }
     810             : INIT:
     811        6748 :   return lfunthetainit_i(data, t, m, bitprec);
     812             : }
     813             : 
     814             : static GEN
     815     4736961 : get_an(GEN an, long n)
     816             : {
     817     4736961 :   if (typ(an) == t_VECSMALL) { long a = an[n]; if (a) return stoi(a); }
     818     4736961 :   else { GEN a = gel(an,n); if (a && !gequal0(a)) return a; }
     819     3421163 :   return NULL;
     820             : }
     821             : /* x * an[n] */
     822             : static GEN
     823    13579080 : mul_an(GEN an, long n, GEN x)
     824             : {
     825    13579080 :   if (typ(an) == t_VECSMALL) { long a = an[n]; if (a) return gmulsg(a,x); }
     826     8500238 :   else { GEN a = gel(an,n); if (a && !gequal0(a)) return gmul(a,x); }
     827     5691320 :   return NULL;
     828             : }
     829             : /* 2*t^a * x **/
     830             : static GEN
     831      151787 : mulT(GEN t, GEN a, GEN x, long prec)
     832             : {
     833      151787 :   if (gequal0(a)) return gmul2n(x,1);
     834       10556 :   return gmul(x, gmul2n(gequal1(a)? t: gpow(t,a,prec), 1));
     835             : }
     836             : 
     837             : static GEN
     838    23388469 : vecan_cmul(void *E, GEN P, long a, GEN x)
     839             : {
     840             :   (void)E;
     841    23388469 :   if (typ(P) == t_VECSMALL)
     842    20433174 :     return (a==0 || !P[a])? NULL: gmulsg(P[a], x);
     843             :   else
     844     2955295 :     return (a==0 || !gel(P,a))? NULL: gmul(gel(P,a), x);
     845             : }
     846             : /* d=2, 2 sum_{n <= N} a(n) (n t)^al q^n, q = exp(-2pi t),
     847             :  * an2[n] = a(n) * n^al */
     848             : static GEN
     849      118090 : theta2_i(GEN an2, long N, GEN t, GEN al, long prec)
     850             : {
     851      118090 :   GEN S, q, pi2 = Pi2n(1,prec);
     852      117930 :   const struct bb_algebra *alg = get_Rg_algebra();
     853      117975 :   setsigne(pi2,-1); q = gexp(gmul(pi2, t), prec);
     854             :   /* Brent-Kung in case the a_n are small integers */
     855      118136 :   S = gen_bkeval(an2, N, q, 1, NULL, alg, vecan_cmul);
     856      118088 :   return mulT(t, al, S, prec);
     857             : }
     858             : static GEN
     859      113763 : theta2(GEN an2, long N, GEN t, GEN al, long prec)
     860             : {
     861      113763 :   pari_sp av = avma;
     862      113763 :   return gerepileupto(av, theta2_i(an2, N, t, al, prec));
     863             : }
     864             : 
     865             : /* d=1, 2 sum_{n <= N} a_n (n t)^al q^(n^2), q = exp(-pi t^2),
     866             :  * an2[n] is a_n n^al */
     867             : static GEN
     868       33698 : theta1(GEN an2, long N, GEN t, GEN al, long prec)
     869             : {
     870       33698 :   GEN q = gexp(gmul(negr(mppi(prec)), gsqr(t)), prec);
     871       33698 :   GEN vexp = gsqrpowers(q, N), S = gen_0;
     872       33698 :   pari_sp av = avma;
     873             :   long n;
     874     6375975 :   for (n = 1; n <= N; n++)
     875             :   {
     876     6342277 :     GEN c = mul_an(an2, n, gel(vexp,n));
     877     6342277 :     if (c)
     878             :     {
     879     5329647 :       S = gadd(S, c);
     880     5329647 :       if (gc_needed(av, 3)) S = gerepileupto(av, S);
     881             :     }
     882             :   }
     883       33698 :   return mulT(t, al, S, prec);
     884             : }
     885             : 
     886             : /* If m > 0, compute m-th derivative of theta(t) = theta0(t/sqrt(N))
     887             :  * with absolute error 2^-bitprec; theta(t)=\sum_{n\ge1}a(n)K(nt/N^(1/2)) */
     888             : GEN
     889       40853 : lfuntheta(GEN data, GEN t, long m, long bitprec)
     890             : {
     891       40853 :   pari_sp ltop = avma;
     892             :   long limt, d;
     893             :   GEN isqN, vecan, Vga, ldata, theta, thetainit, S;
     894       40853 :   long n, prec = nbits2prec(bitprec);
     895       40853 :   t = gprec_w(t, prec);
     896       40853 :   theta = lfunthetacheckinit(data, t, m, bitprec);
     897       40846 :   ldata = linit_get_ldata(theta);
     898       40846 :   thetainit = linit_get_tech(theta);
     899       40846 :   vecan = theta_get_an(thetainit);
     900       40846 :   isqN = theta_get_isqrtN(thetainit);
     901       40846 :   limt = lg(vecan)-1;
     902       40846 :   if (theta == data)
     903       39376 :     limt = minss(limt, lfunthetacost(ldata, t, m, bitprec));
     904       40846 :   if (!limt)
     905             :   {
     906          14 :     set_avma(ltop); S = real_0_bit(-bitprec);
     907          14 :     if (!is_real_t(typ(t)) || !ldata_isreal(ldata))
     908           7 :       S = gerepilecopy(ltop, mkcomplex(S,S));
     909          14 :     return S;
     910             :   }
     911       40832 :   t = gmul(t, isqN);
     912       40832 :   Vga = ldata_get_gammavec(ldata);
     913       40832 :   d = lg(Vga)-1;
     914       40832 :   if (m == 0 && vgaeasytheta(Vga))
     915             :   {
     916       37955 :     if (theta_get_m(thetainit) > 0) vecan = antwist(vecan, Vga, prec);
     917       75910 :     if (d == 1) S = theta1(vecan, limt, t, gel(Vga,1), prec);
     918        4257 :     else        S = theta2_i(vecan, limt, t, vecmin(Vga), prec);
     919             :   }
     920             :   else
     921             :   {
     922        2877 :     GEN K = theta_get_K(thetainit);
     923        2877 :     GEN vroots = mkvroots(d, limt, prec);
     924             :     pari_sp av;
     925        2877 :     t = gpow(t, gdivgs(gen_2,d), prec);
     926        2877 :     S = gen_0; av = avma;
     927     4739838 :     for (n = 1; n <= limt; ++n)
     928             :     {
     929     4736961 :       GEN nt, an = get_an(vecan, n);
     930     4736961 :       if (!an) continue;
     931     1315798 :       nt = gmul(gel(vroots,n), t);
     932     1315798 :       if (m) an = gmul(an, powuu(n, m));
     933     1315798 :       S = gadd(S, gmul(an, gammamellininvrt(K, nt, bitprec)));
     934     1315798 :       if ((n & 0x1ff) == 0) S = gerepileupto(av, S);
     935             :     }
     936        2877 :     if (m) S = gmul(S, gpowgs(isqN, m));
     937             :   }
     938       40832 :   return gerepileupto(ltop, S);
     939             : }
     940             : 
     941             : /*******************************************************************/
     942             : /* Second part: Computation of L-Functions.                        */
     943             : /*******************************************************************/
     944             : 
     945             : struct lfunp {
     946             :   long precmax, Dmax, D, M, m0, nmax, d, vgaell;
     947             :   double k1, dc, dw, dh, MAXs, sub;
     948             :   GEN L, an, bn;
     949             : };
     950             : 
     951             : static void
     952       15552 : lfunparams(GEN ldata, long der, long bitprec, struct lfunp *S)
     953             : {
     954       15552 :   const long derprec = (der > 1)? dbllog2(mpfact(der)): 0; /* log2(der!) */
     955             :   GEN Vga, N, L, k;
     956             :   long k1, d, m, M, flag, nmax;
     957             :   double a, A, E, hd, Ep, d2, suma, maxs, mins, sub, B0,B1;
     958             :   double logN2, logC, Lestimate, Mestimate;
     959             : 
     960       15552 :   Vga = ldata_get_gammavec(ldata);
     961       15552 :   S->d = d = lg(Vga)-1; d2 = d/2.;
     962             : 
     963       15552 :   suma = gtodouble(sumVga(Vga));
     964       15552 :   k = ldata_get_k(ldata);
     965       15552 :   N = ldata_get_conductor(ldata);
     966       15552 :   logN2 = log(gtodouble(N)) / 2;
     967       15552 :   maxs = S->dc + S->dw;
     968       15552 :   mins = S->dc - S->dw;
     969       15552 :   S->MAXs = maxdd(maxs, gtodouble(k)-mins);
     970             : 
     971             :   /* we compute Lambda^(der)(s) / der!; need to compensate for L^(der)(s)
     972             :    * ln |gamma(s)| ~ -(pi/4) \sum_i |Im(s + a_i)|; max with 1: fudge factor */
     973       15552 :   a = (M_PI/(4*M_LN2))*(d*S->dh + sumVgaimpos(Vga));
     974       15552 :   S->D = (long)ceil(bitprec + derprec + maxdd(a, 1));
     975       15552 :   E = M_LN2*S->D; /* D:= required absolute bitprec */
     976             : 
     977       15552 :   Ep = E + maxdd(M_PI * S->dh * d2, (d*S->MAXs + suma - 1) * log(E));
     978       15552 :   hd = d2*M_PI*M_PI / Ep;
     979       15552 :   S->m0 = (long)ceil(M_LN2/hd);
     980       15552 :   hd = M_LN2/S->m0;
     981             : 
     982       15552 :   logC = d2*M_LN2 - log(d2)/2;
     983       15552 :   k1 = ldata_get_k1(ldata);
     984       15552 :   S->k1 = k1; /* assume |a_n| << n^k1 with small implied constant */
     985       15552 :   A = gammavec_expo(d, suma);
     986             : 
     987       15552 :   sub = 0.;
     988       15552 :   if (mins > 1)
     989             :   {
     990        4060 :     GEN sig = dbltor(mins);
     991        4060 :     sub += logN2*mins;
     992        4060 :     if (gammaordinary(Vga, sig))
     993             :     {
     994        3976 :       GEN FVga = gammafactor(Vga);
     995        3976 :       GEN gas = gammafactproduct(FVga, sig, LOWDEFAULTPREC);
     996        3976 :       if (typ(gas) != t_SER)
     997             :       {
     998        3976 :         double dg = dbllog2(gas);
     999        3976 :         if (dg > 0) sub += dg * M_LN2;
    1000             :       }
    1001             :     }
    1002             :   }
    1003       15552 :   S->sub = sub;
    1004       15552 :   M = 1000;
    1005       15552 :   L = cgetg(M+2, t_VECSMALL);
    1006       15552 :   a = S->k1 + A;
    1007             : 
    1008       15552 :   B0 = 5 + E - S->sub + logC + S->k1*logN2; /* 5 extra bits */
    1009       15552 :   B1 = hd * (S->MAXs - S->k1);
    1010       15552 :   Lestimate = dblcoro526(a + S->MAXs - 2./d, d/2.,
    1011       15552 :     E - S->sub + logC - log(2*M_PI*hd) + S->MAXs*logN2);
    1012       15552 :   Mestimate = ((Lestimate > 0? log(Lestimate): 0) + logN2) / hd;
    1013       15552 :   nmax = 0;
    1014       15552 :   flag = 0;
    1015     1544130 :   for (m = 0;; m++)
    1016     1528578 :   {
    1017     1544130 :     double x, H = logN2 - m*hd, B = B0 + m*B1;
    1018             :     long n;
    1019     1544130 :     x = dblcoro526(a, d/2., B);
    1020     1544130 :     n = floor(x*exp(H));
    1021     1544130 :     if (n > nmax) nmax = n;
    1022     1544130 :     if (m > M) { M *= 2; L = vecsmall_lengthen(L,M+2); }
    1023     1544130 :     L[m+1] = n;
    1024     1544130 :     if (n == 0) { if (++flag > 2 && m > Mestimate) break; } else flag = 0;
    1025             :   }
    1026       15552 :   m -= 2; while (m > 0 && !L[m]) m--;
    1027       15552 :   if (m == 0) { nmax = 1; L[1] = 1; m = 1; } /* can happen for tiny bitprec */
    1028       15552 :   setlg(L, m+1); S->M = m-1;
    1029       15552 :   S->L = L;
    1030       15552 :   S->nmax = nmax;
    1031             : 
    1032       15552 :   S->Dmax = S->D + (long)ceil((S->M * hd * S->MAXs - S->sub) / M_LN2);
    1033       15552 :   if (S->Dmax < S->D) S->Dmax = S->D;
    1034       15552 :   S->precmax = nbits2prec(S->Dmax);
    1035       15552 :   if (DEBUGLEVEL > 1)
    1036           0 :     err_printf("Dmax=%ld, D=%ld, M = %ld, nmax = %ld, m0 = %ld\n",
    1037             :                S->Dmax,S->D,S->M,S->nmax, S->m0);
    1038       15552 : }
    1039             : 
    1040             : static GEN
    1041        4662 : lfuninit_pol(GEN v, GEN poqk, long prec)
    1042             : {
    1043        4662 :   long m, M = lg(v) - 2;
    1044        4662 :   GEN pol = cgetg(M+3, t_POL);
    1045        4662 :   pol[1] = evalsigne(1) | evalvarn(0);
    1046        4662 :   gel(pol, 2) = gprec_w(gmul2n(gel(v,1), -1), prec);
    1047      281974 :   for (m = 2; m <= M+1; m++)
    1048      277312 :     gel(pol, m+1) = gprec_w(gmul(gel(poqk,m), gel(v,m)), prec);
    1049        4662 :   return RgX_renormalize_lg(pol, M+3);
    1050             : }
    1051             : 
    1052             : static void
    1053       50790 : worker_init(long q, GEN *an, GEN *bn, GEN *AB, GEN *A, GEN *B)
    1054             : {
    1055       50790 :   if (typ(*bn) == t_INT) *bn = NULL;
    1056       50790 :   if (*bn)
    1057             :   {
    1058        1015 :     *AB = cgetg(3, t_VEC);
    1059        1014 :     gel(*AB,1) = *A = cgetg(q+1, t_VEC);
    1060        1014 :     gel(*AB,2) = *B = cgetg(q+1, t_VEC);
    1061        1014 :     if (typ(an) == t_VEC) *an = RgV_kill0(*an);
    1062        1014 :     if (typ(bn) == t_VEC) *bn = RgV_kill0(*bn);
    1063             :   }
    1064             :   else
    1065             :   {
    1066       49775 :     *B = NULL;
    1067       49775 :     *AB = *A = cgetg(q+1, t_VEC);
    1068       49748 :     if (typ(*an) == t_VEC) *an = RgV_kill0(*an);
    1069             :   }
    1070       50386 : }
    1071             : GEN
    1072       10260 : lfuninit_theta2_worker(long r, GEN L, GEN qk, GEN a, GEN di, GEN an, GEN bn)
    1073             : {
    1074       10260 :   long q, m, prec = di[1], M = di[2], m0 = di[3], L0 = lg(an)-1;
    1075             :   GEN AB, A, B;
    1076       10260 :   worker_init((M - r) / m0 + 1, &an, &bn, &AB, &A, &B);
    1077      116838 :   for (q = 0, m = r; m <= M; m += m0, q++)
    1078             :   {
    1079      106708 :     GEN t = gel(qk, m+1);
    1080      106708 :     long N = minss(L[m+1],L0);
    1081      106764 :     gel(A, q+1) = theta2(an, N, t, a, prec); /* theta(exp(mh)) */
    1082      106855 :     if (bn) gel(B, q+1) = theta2(bn, N, t, a, prec);
    1083             :   }
    1084       10130 :   return AB;
    1085             : }
    1086             : 
    1087             : /* theta(exp(mh)) ~ sum_{n <= N} a(n) k[m,n] */
    1088             : static GEN
    1089      168518 : an_msum(GEN an, long N, GEN vKm)
    1090             : {
    1091      168518 :   pari_sp av = avma;
    1092      168518 :   GEN s = gen_0;
    1093             :   long n;
    1094     7402255 :   for (n = 1; n <= N; n++)
    1095             :   {
    1096     7234813 :     GEN c = mul_an(an, n, gel(vKm,n));
    1097     7234976 :     if (c) s = gadd(s, c);
    1098             :   }
    1099      167442 :   return gerepileupto(av, s);
    1100             : }
    1101             : 
    1102             : GEN
    1103       40523 : lfuninit_worker(long r, GEN K, GEN L, GEN peh2d, GEN vroots, GEN dr, GEN di,
    1104             :                 GEN an, GEN bn)
    1105             : {
    1106       40523 :   pari_sp av0 = avma;
    1107       40523 :   long m, n, q, L0 = lg(an)-1;
    1108       40523 :   double sig0 = rtodbl(gel(dr,1));
    1109       40545 :   double sub2 = rtodbl(gel(dr,1));
    1110       40531 :   double k1 = rtodbl(gel(dr,3));
    1111       40528 :   double MAXs = rtodbl(gel(dr,4));
    1112       40519 :   long D = di[1], M = di[2], m0 = di[3];
    1113       40519 :   double M0 = sig0? sub2 / sig0: 1./0.;
    1114       40519 :   GEN AB, A, B, vK = cgetg(M/m0 + 2, t_VEC);
    1115             : 
    1116      206718 :   for (q = 0, m = r; m <= M; m += m0, q++)
    1117      166232 :     gel(vK, q+1) = const_vec(L[m+1], NULL);
    1118       40486 :   worker_init(q, &an, &bn, &AB, &A, &B);
    1119      207181 :   for (m -= m0, q--; m >= 0; m -= m0, q--)
    1120             :   {
    1121      166636 :     double c1 = D + ((m > M0)? m * sig0 - sub2 : 0);
    1122      166636 :     GEN vKm = gel(vK,q+1); /* conceptually K(m,n) */
    1123     7410068 :     for (n = 1; n <= L[m+1]; n++)
    1124             :     {
    1125             :       GEN t2d, kmn;
    1126     7243366 :       long nn, mm, qq, p = 0;
    1127             :       double c, c2;
    1128             :       pari_sp av;
    1129             : 
    1130     7243366 :       if (gel(vKm, n)) continue; /* done already */
    1131     5390421 :       c = c1 + k1 * log2(n);
    1132             :       /* n *= 2; m -= m0 => c += c2, provided m >= M0. Else c += k1 */
    1133     5390421 :       c2 = k1 - MAXs;
    1134             :       /* p = largest (absolute) accuracy to which we need K(m,n) */
    1135    15422030 :       for (mm=m,nn=n; mm >= M0;)
    1136             :       {
    1137     9348994 :         if (nn <= L[mm+1] && (gel(an, nn) || (bn && gel(bn, nn))))
    1138     2306872 :           p = maxuu(p, (ulong)c);
    1139     9350006 :         nn <<= 1;
    1140     9350006 :         mm -= m0; if (mm >= M0) c += c2; else { c += k1; break; }
    1141             :       }
    1142             :       /* mm < M0 || nn > L[mm+1] */
    1143     6767619 :       for (         ; mm >= 0; nn<<=1,mm-=m0,c+=k1)
    1144     1376186 :         if (nn <= L[mm+1] && (gel(an, nn) || (bn && gel(bn, nn))))
    1145      256421 :           p = maxuu(p, (ulong)c);
    1146     5391433 :       if (!p) continue; /* a_{n 2^v} = 0 for all v in range */
    1147     1829020 :       av = avma;
    1148     1829020 :       t2d = mpmul(gel(vroots,n), gel(peh2d,m+1));/*(n exp(mh)/sqrt(N))^(2/d)*/
    1149     1828968 :       kmn = gerepileupto(av, gammamellininvrt(K, t2d, p));
    1150     5630166 :       for (qq=q,mm=m,nn=n; mm >= 0; nn<<=1,mm-=m0,qq--)
    1151     3802092 :         if (nn <= L[mm+1]) gmael(vK, qq+1, nn) = kmn;
    1152             :     }
    1153             :   }
    1154      207085 :   for (q = 0, m = r; m <= M; m += m0, q++)
    1155             :   {
    1156      166569 :     long N = minss(L0, L[m+1]);
    1157      166561 :     gel(A, q+1) = an_msum(an, N, gel(vK,q+1));
    1158      166539 :     if (bn) gel(B, q+1) = an_msum(bn, N, gel(vK,q+1));
    1159             :   }
    1160       40516 :   return gerepileupto(av0, AB);
    1161             : }
    1162             : /* return A = [\theta(exp(mh)), m=0..M], theta(t) = sum a(n) K(n/sqrt(N) t),
    1163             :  * h = log(2)/m0. If bn != NULL, return the pair [A, B] */
    1164             : static GEN
    1165        4508 : lfuninit_ab(GEN theta, GEN h, struct lfunp *S)
    1166             : {
    1167        4508 :   const long M = S->M, prec = S->precmax;
    1168        4508 :   GEN tech = linit_get_tech(theta), isqN = theta_get_isqrtN(tech);
    1169        4508 :   GEN an = S->an, bn = S->bn, va, vb;
    1170             :   struct pari_mt pt;
    1171             :   GEN worker;
    1172             :   long m0, r, pending;
    1173             : 
    1174        4508 :   if (S->vgaell)
    1175             :   { /* d=2 and Vga = [a,a+1] */
    1176        1981 :     GEN a = vecmin(ldata_get_gammavec(linit_get_ldata(theta)));
    1177        1981 :     GEN qk = gpowers0(mpexp(h), M, isqN);
    1178        1981 :     m0 = minss(M+1, pari_mt_nbthreads);
    1179        1981 :     worker = snm_closure(is_entry("_lfuninit_theta2_worker"),
    1180             :                          mkvecn(6, S->L, qk, a, mkvecsmall3(prec, M, m0),
    1181             :                                 an, bn? bn: gen_0));
    1182             :   }
    1183             :   else
    1184             :   {
    1185             :     GEN vroots, peh2d, d2;
    1186        2527 :     double sig0 = S->MAXs / S->m0, sub2 = S->sub / M_LN2;
    1187             :     /* For all 0<= m <= M, and all n <= L[m+1] such that a_n!=0, we must compute
    1188             :      *   k[m,n] = K(n exp(mh)/sqrt(N))
    1189             :      * with ln(absolute error) <= E + max(mh sigma - sub, 0) + k1 * log(n).
    1190             :      * N.B. we use the 'rt' variant and pass argument (n exp(mh)/sqrt(N))^(2/d).
    1191             :      * Speedup: if n' = 2n and m' = m - m0 >= 0; then k[m,n] = k[m',n']. */
    1192             :     /* vroots[n] = n^(2/d) */
    1193        2527 :     vroots = mkvroots(S->d, S->nmax, prec);
    1194        2527 :     d2 = gdivgs(gen_2, S->d);
    1195             :     /* peh2d[m+1] = (exp(mh)/sqrt(N))^(2/d) */
    1196        2527 :     peh2d = gpowers0(gexp(gmul(d2,h), prec), M, gpow(isqN, d2, prec));
    1197        2527 :     m0 = S->m0;
    1198        2527 :     worker = snm_closure(is_entry("_lfuninit_worker"),
    1199             :                          mkvecn(8, theta_get_K(tech), S->L, peh2d, vroots,
    1200             :                                 mkvec4(dbltor(sig0), dbltor(sub2),
    1201             :                                        dbltor(S->k1), dbltor(S->MAXs)),
    1202             :                                 mkvecsmall3(S->D, M, m0),
    1203             :                                 an, bn? bn: gen_0));
    1204             :     /* For each 0 <= m <= M, we will sum for n<=L[m+1] a(n) K(m,n)
    1205             :      * bit accuracy for K(m,n): D + k1*log2(n) + 1_{m > M0} (m*sig0 - sub2)
    1206             :      * We restrict m to arithmetic progressions r mod m0 to save memory and
    1207             :      * allow parallelization */
    1208             :   }
    1209        4508 :   va = cgetg(M+2, t_VEC);
    1210        4508 :   vb = bn? cgetg(M+2, t_VEC): NULL;
    1211        4508 :   mt_queue_start_lim(&pt, worker, m0);
    1212        4508 :   pending = 0;
    1213       68076 :   for (r = 0; r < m0 || pending; r++)
    1214             :   { /* m = q m0 + r */
    1215             :     GEN done, A, B;
    1216             :     long q, m, workid;
    1217       63568 :     mt_queue_submit(&pt, r, r < m0 ? mkvec(utoi(r)): NULL);
    1218       63568 :     done = mt_queue_get(&pt, &workid, &pending);
    1219       63568 :     if (!done) continue;
    1220       50915 :     if (bn) { A = gel(done,1); B = gel(done,2); } else { A = done; B = NULL; }
    1221      324909 :     for (q = 0, m = workid; m <= M; m += m0, q++)
    1222             :     {
    1223      273994 :       gel(va, m+1) = gel(A, q+1);
    1224      273994 :       if (bn) gel(vb, m+1) = gel(B, q+1);
    1225             :     }
    1226             :   }
    1227        4508 :   mt_queue_end(&pt);
    1228        4508 :   return bn? mkvec2(va, vb): va;
    1229             : }
    1230             : 
    1231             : static void
    1232       81578 : parse_dom(double k, GEN dom, struct lfunp *S)
    1233             : {
    1234       81578 :   long l = lg(dom);
    1235       81578 :   if (typ(dom)!=t_VEC) pari_err_TYPE("lfuninit [domain]", dom);
    1236       81578 :   if (l == 2)
    1237             :   {
    1238       44002 :     S->dc = k/2.;
    1239       44002 :     S->dw = 0.;
    1240       44002 :     S->dh = gtodouble(gel(dom,1));
    1241             :   }
    1242       37576 :   else if (l == 3)
    1243             :   {
    1244         301 :     S->dc = k/2.;
    1245         301 :     S->dw = gtodouble(gel(dom,1));
    1246         301 :     S->dh = gtodouble(gel(dom,2));
    1247             :   }
    1248       37275 :   else if (l == 4)
    1249             :   {
    1250       37275 :     S->dc = gtodouble(gel(dom,1));
    1251       37275 :     S->dw = gtodouble(gel(dom,2));
    1252       37275 :     S->dh = gtodouble(gel(dom,3));
    1253             :   }
    1254             :   else
    1255             :   {
    1256           0 :     pari_err_TYPE("lfuninit [domain]", dom);
    1257           0 :     S->dc = S->dw = S->dh = 0; /*-Wall*/
    1258             :   }
    1259       81578 :   if (S->dw < 0 || S->dh < 0) pari_err_TYPE("lfuninit [domain]", dom);
    1260       81578 : }
    1261             : 
    1262             : /* do we have dom \subset dom0 ? dom = [center, width, height] */
    1263             : int
    1264       14110 : sdomain_isincl(double k, GEN dom, GEN dom0)
    1265             : {
    1266             :   struct lfunp S0, S;
    1267       14110 :   parse_dom(k, dom, &S);
    1268       14110 :   parse_dom(k, dom0, &S0);
    1269       14110 :   return S0.dc - S0.dw <= S.dc - S.dw
    1270       14110 :       && S0.dc + S0.dw >= S.dc + S.dw && S0.dh >= S.dh;
    1271             : }
    1272             : 
    1273             : static int
    1274       14110 : checklfuninit(GEN linit, GEN dom, long der, long bitprec)
    1275             : {
    1276       14110 :   GEN ldata = linit_get_ldata(linit);
    1277       14110 :   GEN domain = lfun_get_domain(linit_get_tech(linit));
    1278       14110 :   return domain_get_der(domain) >= der
    1279       14110 :     && domain_get_bitprec(domain) >= bitprec
    1280       28220 :     && sdomain_isincl(gtodouble(ldata_get_k(ldata)), dom, domain_get_dom(domain));
    1281             : }
    1282             : 
    1283             : GEN
    1284        5194 : lfuninit_make(long t, GEN ldata, GEN molin, GEN domain)
    1285             : {
    1286        5194 :   GEN Vga = ldata_get_gammavec(ldata);
    1287        5194 :   long d = lg(Vga)-1;
    1288        5194 :   GEN hardy, w2 = gen_1, k2 = gmul2n(ldata_get_k(ldata), -1);
    1289        5194 :   GEN expot = gdivgs(gadd(gmulsg(d, gsubgs(k2, 1)), sumVga(Vga)), 4);
    1290        5194 :   if (typ(ldata_get_dual(ldata))==t_INT)
    1291             :   {
    1292        5040 :     GEN eno = ldata_get_rootno(ldata);
    1293        5040 :     long prec = nbits2prec( domain_get_bitprec(domain) );
    1294        5040 :     if (!isint1(eno)) w2 = ginv(gsqrt(eno, prec));
    1295             :   }
    1296        5194 :   hardy = mkvec4(k2, w2, expot, gammafactor(Vga));
    1297        5194 :   return mkvec3(mkvecsmall(t),ldata, mkvec3(domain, molin, hardy));
    1298             : }
    1299             : 
    1300             : static void
    1301        2527 : lfunparams2(struct lfunp *S)
    1302             : {
    1303        2527 :   GEN L = S->L, an = S->an, bn = S->bn;
    1304             :   double pmax;
    1305        2527 :   long m, nan, nmax, neval, M = S->M;
    1306             : 
    1307        2527 :   S->vgaell = 0;
    1308             :   /* try to reduce parameters now we know the a_n (some may be 0) */
    1309        2527 :   if (typ(an) == t_VEC) an = RgV_kill0(an);
    1310        2527 :   nan = S->nmax; /* lg(an)-1 may be large than this */
    1311        2527 :   nmax = neval = 0;
    1312        2527 :   if (!bn)
    1313      168315 :     for (m = 0; m <= M; m++)
    1314             :     {
    1315      165809 :       long n = minss(nan, L[m+1]);
    1316      165809 :       while (n > 0 && !gel(an,n)) n--;
    1317      165809 :       if (n > nmax) nmax = n;
    1318      165809 :       neval += n;
    1319      165809 :       L[m+1] = n; /* reduce S->L[m+1] */
    1320             :     }
    1321             :   else
    1322             :   {
    1323          21 :     if (typ(bn) == t_VEC) bn = RgV_kill0(bn);
    1324        1036 :     for (m = 0; m <= M; m++)
    1325             :     {
    1326        1015 :       long n = minss(nan, L[m+1]);
    1327        1015 :       while (n > 0 && !gel(an,n) && !gel(bn,n)) n--;
    1328        1015 :       if (n > nmax) nmax = n;
    1329        1015 :       neval += n;
    1330        1015 :       L[m+1] = n; /* reduce S->L[m+1] */
    1331             :     }
    1332             :   }
    1333        2527 :   if (DEBUGLEVEL >= 1) err_printf("expected evaluations: %ld\n", neval);
    1334        2527 :   for (; M > 0; M--)
    1335        2527 :     if (L[M+1]) break;
    1336        2527 :   setlg(L, M+2);
    1337        2527 :   S->M = M;
    1338        2527 :   S->nmax = nmax;
    1339             : 
    1340             :   /* need K(n*exp(mh)/sqrt(N)) to absolute accuracy
    1341             :    *   D + k1*log(n) + max(m * sig0 - sub2, 0) */
    1342        2527 :   pmax = S->D + S->k1 * log2(L[1]);
    1343        2527 :   if (S->MAXs)
    1344             :   {
    1345        2527 :     double sig0 = S->MAXs/S->m0, sub2 = S->sub / M_LN2;
    1346      120687 :     for (m = ceil(sub2 / sig0); m <= S->M; m++)
    1347             :     {
    1348      118160 :       double c = S->D + m*sig0 - sub2;
    1349      118160 :       if (S->k1 > 0) c += S->k1 * log2(L[m+1]);
    1350      118160 :       pmax = maxdd(pmax, c);
    1351             :     }
    1352             :   }
    1353        2527 :   S->Dmax = pmax;
    1354        2527 :   S->precmax = nbits2prec(pmax);
    1355        2527 : }
    1356             : 
    1357             : static GEN
    1358        4515 : lfun_init_theta(GEN ldata, GEN eno, struct lfunp *S)
    1359             : {
    1360        4515 :   GEN an2, dual, tdom = NULL, Vga = ldata_get_gammavec(ldata);
    1361        4515 :   long L, prec = S->precmax;
    1362        4515 :   if (eno)
    1363        3430 :     L = S->nmax;
    1364             :   else
    1365             :   {
    1366        1085 :     tdom = dbltor(sqrt(0.5));
    1367        1085 :     L = maxss(S->nmax, lfunthetacost(ldata, tdom, 0, S->D));
    1368             :   }
    1369        4515 :   dual = ldata_get_dual(ldata);
    1370        4515 :   S->an = ldata_vecan(ldata_get_an(ldata), L, prec);
    1371        4508 :   S->bn = typ(dual)==t_INT? NULL: ldata_vecan(dual, S->nmax, prec);
    1372        4508 :   if (!vgaell(Vga)) lfunparams2(S);
    1373             :   else
    1374             :   {
    1375        1981 :     S->an = antwist(S->an, Vga, prec);
    1376        1981 :     if (S->bn) S->bn = antwist(S->bn, Vga, prec);
    1377        1981 :     S->vgaell = 1;
    1378             :   }
    1379        4508 :   an2 = lg(Vga)-1 == 1? antwist(S->an, Vga, prec): S->an;
    1380        4508 :   return lfunthetainit0(ldata, tdom, an2, 0, S->Dmax, 0);
    1381             : }
    1382             : 
    1383             : GEN
    1384       11037 : lfuncost(GEN L, GEN dom, long der, long bitprec)
    1385             : {
    1386       11037 :   pari_sp av = avma;
    1387       11037 :   GEN ldata = lfunmisc_to_ldata_shallow(L);
    1388       11037 :   GEN k = ldata_get_k(ldata);
    1389             :   struct lfunp S;
    1390             : 
    1391       11037 :   parse_dom(gtodouble(k), dom, &S);
    1392       11037 :   lfunparams(ldata, der, bitprec, &S);
    1393       11037 :   set_avma(av); return mkvecsmall2(S.nmax, S.Dmax);
    1394             : }
    1395             : GEN
    1396          42 : lfuncost0(GEN L, GEN dom, long der, long bitprec)
    1397             : {
    1398          42 :   pari_sp av = avma;
    1399             :   GEN C;
    1400             : 
    1401          42 :   if (is_linit(L))
    1402             :   {
    1403          28 :     GEN tech = linit_get_tech(L);
    1404          28 :     GEN domain = lfun_get_domain(tech);
    1405          28 :     dom = domain_get_dom(domain);
    1406          28 :     der = domain_get_der(domain);
    1407          28 :     bitprec = domain_get_bitprec(domain);
    1408          28 :     if (linit_get_type(L) == t_LDESC_PRODUCT)
    1409             :     {
    1410          21 :       GEN v = lfunprod_get_fact(linit_get_tech(L)), F = gel(v,1);
    1411          21 :       long i, l = lg(F);
    1412          21 :       C = cgetg(l, t_VEC);
    1413          77 :       for (i = 1; i < l; ++i)
    1414          56 :         gel(C, i) = zv_to_ZV( lfuncost(gel(F,i), dom, der, bitprec) );
    1415          21 :       return gerepileupto(av, C);
    1416             :     }
    1417             :   }
    1418          21 :   if (!dom) pari_err_TYPE("lfuncost [missing s domain]", L);
    1419          21 :   C = lfuncost(L,dom,der,bitprec);
    1420          21 :   return gerepileupto(av, zv_to_ZV(C));
    1421             : }
    1422             : 
    1423             : GEN
    1424       19101 : lfuninit(GEN lmisc, GEN dom, long der, long bitprec)
    1425             : {
    1426       19101 :   pari_sp ltop = avma;
    1427             :   GEN poqk, AB, R, h, theta, ldata, eno, r, domain, molin, k;
    1428             :   struct lfunp S;
    1429             : 
    1430       19101 :   if (is_linit(lmisc))
    1431             :   {
    1432       14159 :     long t = linit_get_type(lmisc);
    1433       14159 :     if (t==t_LDESC_INIT || t==t_LDESC_PRODUCT)
    1434             :     {
    1435       14110 :       if (checklfuninit(lmisc, dom, der, bitprec)) return lmisc;
    1436          21 :       pari_warn(warner,"lfuninit: insufficient initialization");
    1437             :     }
    1438             :   }
    1439        5012 :   ldata = lfunmisc_to_ldata_shallow(lmisc);
    1440             : 
    1441        5012 :   if (ldata_get_type(ldata)==t_LFUN_NF)
    1442             :   {
    1443         497 :     GEN T = gel(ldata_get_an(ldata), 2);
    1444         497 :     return lfunzetakinit(T, dom, der, 0, bitprec);
    1445             :   }
    1446        4515 :   k = ldata_get_k(ldata);
    1447        4515 :   parse_dom(gtodouble(k), dom, &S);
    1448        4515 :   lfunparams(ldata, der, bitprec, &S);
    1449        4515 :   r = ldata_get_residue(ldata);
    1450             :   /* Note: all guesses should already have been performed (thetainit more
    1451             :    * expensive than needed: should be either tdom = 1 or bitprec = S.D).
    1452             :    * BUT if the root number / polar part do not have an algebraic
    1453             :    * expression, there is no way to do this until we know the
    1454             :    * precision, i.e. now. So we can't remove guessing code from here and
    1455             :    * lfun_init_theta */
    1456        4515 :   if (r && isintzero(r)) eno = NULL;
    1457             :   else
    1458             :   {
    1459        4515 :     eno = ldata_get_rootno(ldata);
    1460        4515 :     if (isintzero(eno)) eno = NULL;
    1461             :   }
    1462        4515 :   theta = lfun_init_theta(ldata, eno, &S);
    1463        4508 :   if (eno && lg(ldata)==7)
    1464        2037 :     R = gen_0;
    1465             :   else
    1466             :   {
    1467        2471 :     GEN v = lfunrootres(theta, S.D);
    1468        2471 :     ldata = shallowcopy(ldata);
    1469        2471 :     gel(ldata, 6) = gel(v,3);
    1470        2471 :     r = gel(v,1);
    1471        2471 :     if (isintzero(r))
    1472        1071 :       setlg(ldata,7); /* no pole */
    1473             :     else
    1474        1400 :       gel(ldata, 7) = r;
    1475        2471 :     R = lfunrtoR(ldata, nbits2prec(S.D));
    1476             :   }
    1477        4508 :   h = divru(mplog2(S.precmax), S.m0);
    1478             :   /* exp(kh/2 . [0..M]) */
    1479        4508 :   poqk = gpowers(gprec_w(mpexp(gmul2n(gmul(k,h), -1)), S.precmax), S.M);
    1480        4508 :   AB = lfuninit_ab(theta, h, &S);
    1481        4508 :   if (S.bn)
    1482             :   {
    1483         154 :     GEN A = gel(AB,1), B = gel(AB,2);
    1484         154 :     A = lfuninit_pol(A, poqk, S.precmax);
    1485         154 :     B = lfuninit_pol(B, poqk, S.precmax);
    1486         154 :     AB = mkvec2(A, B);
    1487             :   }
    1488             :   else
    1489        4354 :     AB = lfuninit_pol(AB, poqk, S.precmax);
    1490        4508 :   molin = mkvec3(h, AB, R);
    1491        4508 :   domain = mkvec2(dom, mkvecsmall2(der, bitprec));
    1492        4508 :   return gerepilecopy(ltop, lfuninit_make(t_LDESC_INIT, ldata, molin, domain));
    1493             : }
    1494             : 
    1495             : GEN
    1496         406 : lfuninit0(GEN lmisc, GEN dom, long der, long bitprec)
    1497             : {
    1498         406 :   GEN z = lfuninit(lmisc, dom, der, bitprec);
    1499         406 :   return z == lmisc? gcopy(z): z;
    1500             : }
    1501             : 
    1502             : /* If s is a pole of Lambda, return polar part at s; else return NULL */
    1503             : static GEN
    1504        4247 : lfunpoleresidue(GEN R, GEN s)
    1505             : {
    1506             :   long j;
    1507       11873 :   for (j = 1; j < lg(R); j++)
    1508             :   {
    1509        8158 :     GEN Rj = gel(R, j), be = gel(Rj, 1);
    1510        8158 :     if (gequal(s, be)) return gel(Rj, 2);
    1511             :   }
    1512        3715 :   return NULL;
    1513             : }
    1514             : 
    1515             : /* Compute contribution of polar part at s when not a pole. */
    1516             : static GEN
    1517        6681 : veccothderivn(GEN a, long n)
    1518             : {
    1519             :   long i;
    1520        6681 :   pari_sp av = avma;
    1521        6681 :   GEN c = pol_x(0), cp = mkpoln(3, gen_m1, gen_0, gen_1);
    1522        6681 :   GEN v = cgetg(n+2, t_VEC);
    1523        6681 :   gel(v, 1) = poleval(c, a);
    1524       20106 :   for(i = 2; i <= n+1; i++)
    1525             :   {
    1526       13425 :     c = ZX_mul(ZX_deriv(c), cp);
    1527       13425 :     gel(v, i) = gdiv(poleval(c, a), mpfact(i-1));
    1528             :   }
    1529        6681 :   return gerepilecopy(av, v);
    1530             : }
    1531             : 
    1532             : static GEN
    1533        6744 : polepart(long n, GEN h, GEN C)
    1534             : {
    1535        6744 :   GEN h2n = gpowgs(gdiv(h, gen_2), n-1);
    1536        6744 :   GEN res = gmul(h2n, gel(C,n));
    1537        6744 :   return odd(n)? res : gneg(res);
    1538             : }
    1539             : 
    1540             : static GEN
    1541        3337 : lfunsumcoth(GEN R, GEN s, GEN h, long prec)
    1542             : {
    1543             :   long i,j;
    1544        3337 :   GEN S = gen_0;
    1545       10018 :   for (j = 1; j < lg(R); ++j)
    1546             :   {
    1547        6681 :     GEN r = gel(R,j), be = gel(r,1), Rj = gel(r, 2);
    1548        6681 :     long e = valp(Rj);
    1549        6681 :     GEN z1 = gexpm1(gmul(h, gsub(s,be)), prec); /* exp(h(s-beta))-1 */
    1550        6681 :     GEN c1 = gaddgs(gdivsg(2, z1), 1); /* coth((h/2)(s-beta)) */
    1551        6681 :     GEN C1 = veccothderivn(c1, 1-e);
    1552       13425 :     for (i = e; i < 0; i++)
    1553             :     {
    1554        6744 :       GEN Rbe = mysercoeff(Rj, i);
    1555        6744 :       GEN p1 = polepart(-i, h, C1);
    1556        6744 :       S = gadd(S, gmul(Rbe, p1));
    1557             :     }
    1558             :   }
    1559        3337 :   return gmul2n(S, -1);
    1560             : }
    1561             : 
    1562             : static GEN lfunlambda_OK(GEN linit, GEN s, GEN sdom, long bitprec);
    1563             : /* L is a t_LDESC_PRODUCT Linit */
    1564             : static GEN
    1565        1279 : lfunlambda_product(GEN L, GEN s, GEN sdom, long bitprec)
    1566             : {
    1567        1279 :   GEN ldata = linit_get_ldata(L), v = lfunprod_get_fact(linit_get_tech(L));
    1568        1279 :   GEN r = gen_1, F = gel(v,1), E = gel(v,2), C = gel(v,3), cs = conj_i(s);
    1569        1279 :   long i, l = lg(F), isreal = gequal(imag_i(s), imag_i(cs));
    1570        4439 :   for (i = 1; i < l; ++i)
    1571             :   {
    1572        3160 :     GEN f = lfunlambda_OK(gel(F, i), s, sdom, bitprec);
    1573        3160 :     if (E[i])
    1574        3160 :       r = gmul(r, gpowgs(f, E[i]));
    1575        3160 :     if (C[i])
    1576             :     {
    1577         378 :       GEN fc = isreal? f: lfunlambda_OK(gel(F, i), cs, sdom, bitprec);
    1578         378 :       r = gmul(r, gpowgs(conj_i(fc), C[i]));
    1579             :     }
    1580             :   }
    1581        1279 :   return (ldata_isreal(ldata) && gequal0(imag_i(s)))? real_i(r): r;
    1582             : }
    1583             : 
    1584             : /* s a t_SER */
    1585             : static long
    1586        1107 : der_level(GEN s)
    1587        1107 : { return signe(s)? lg(s)-3: valp(s)-1; }
    1588             : 
    1589             : /* s a t_SER; return coeff(s, X^0) */
    1590             : static GEN
    1591         217 : ser_coeff0(GEN s) { return simplify_shallow(polcoef_i(s, 0, -1)); }
    1592             : 
    1593             : static GEN
    1594        5404 : get_domain(GEN s, GEN *dom, long *der)
    1595             : {
    1596        5404 :   GEN sa = s;
    1597        5404 :   *der = 0;
    1598        5404 :   switch(typ(s))
    1599             :   {
    1600             :     case t_POL:
    1601           7 :     case t_RFRAC: s = toser_i(s);
    1602             :     case t_SER:
    1603         217 :       *der = der_level(s);
    1604         217 :       sa = ser_coeff0(s);
    1605             :   }
    1606        5404 :   *dom = mkvec3(real_i(sa), gen_0, gabs(imag_i(sa),DEFAULTPREC));
    1607        5404 :   return s;
    1608             : }
    1609             : 
    1610             : /* assume lmisc is an linit, s went through get_domain and s/bitprec belong
    1611             :  * to domain */
    1612             : static GEN
    1613       20182 : lfunlambda_OK(GEN linit, GEN s, GEN sdom, long bitprec)
    1614             : {
    1615             :   GEN eno, ldata, tech, h, pol;
    1616       20182 :   GEN S, S0 = NULL, k2, cost;
    1617             :   long prec, prec0;
    1618             :   struct lfunp D, D0;
    1619             : 
    1620       20182 :   if (linit_get_type(linit) == t_LDESC_PRODUCT)
    1621        1279 :     return lfunlambda_product(linit, s, sdom, bitprec);
    1622       18903 :   ldata = linit_get_ldata(linit);
    1623       18903 :   eno = ldata_get_rootno(ldata);
    1624       18903 :   tech = linit_get_tech(linit);
    1625       18903 :   h = lfun_get_step(tech); prec = realprec(h);
    1626             :   /* try to reduce accuracy */
    1627       18903 :   parse_dom(0, sdom, &D0);
    1628       18903 :   parse_dom(0, domain_get_dom(lfun_get_domain(tech)), &D);
    1629       18903 :   if (0.8 * D.dh > D0.dh)
    1630             :   {
    1631       10960 :     cost = lfuncost(linit, sdom, typ(s)==t_SER? der_level(s): 0, bitprec);
    1632       10960 :     prec0 = nbits2prec(cost[2]);
    1633       10960 :     if (prec0 < prec) { prec = prec0; h = gprec_w(h, prec); }
    1634             :   }
    1635       18903 :   pol = lfun_get_pol(tech);
    1636       18903 :   s = gprec_w(s, prec);
    1637       18903 :   if (ldata_get_residue(ldata))
    1638             :   {
    1639        3778 :     GEN R = lfun_get_Residue(tech);
    1640        3778 :     GEN Ra = lfunpoleresidue(R, s);
    1641        3778 :     if (Ra) return gprec_w(Ra, nbits2prec(bitprec));
    1642        3337 :     S0 = lfunsumcoth(R, s, h, prec);
    1643             :   }
    1644       18462 :   k2 = lfun_get_k2(tech);
    1645       18462 :   if (typ(pol)==t_POL && typ(s) != t_SER && gequal(real_i(s), k2))
    1646       13065 :   { /* on critical line: shortcut */
    1647       13065 :     GEN polz, b = imag_i(s);
    1648       13065 :     polz = gequal0(b)? poleval(pol,gen_1): poleval(pol, expIr(gmul(h,b)));
    1649       13065 :     S = gadd(polz, gmul(eno, conj_i(polz)));
    1650             :   }
    1651             :   else
    1652             :   {
    1653        5397 :     GEN z = gexp(gmul(h, gsub(s, k2)), prec);
    1654        5397 :     GEN zi = ginv(z), zc = conj_i(zi);
    1655        5397 :     if (typ(pol)==t_POL)
    1656        5201 :       S = gadd(poleval(pol, z), gmul(eno, conj_i(poleval(pol, zc))));
    1657             :     else
    1658         196 :       S = gadd(poleval(gel(pol,1), z), gmul(eno, poleval(gel(pol,2), zi)));
    1659             :   }
    1660       18462 :   if (S0) S = gadd(S,S0);
    1661       18462 :   return gprec_w(gmul(S,h), nbits2prec(bitprec));
    1662             : }
    1663             : GEN
    1664         882 : lfunlambda(GEN lmisc, GEN s, long bitprec)
    1665             : {
    1666         882 :   pari_sp av = avma;
    1667             :   GEN linit, dom, z;
    1668             :   long der;
    1669         882 :   s = get_domain(s, &dom, &der);
    1670         882 :   linit = lfuninit(lmisc, dom, der, bitprec);
    1671         882 :   z = lfunlambda_OK(linit,s, dom, bitprec);
    1672         882 :   return gerepilecopy(av, z);
    1673             : }
    1674             : 
    1675             : /* assume lmisc is an linit, s went through get_domain and s/bitprec belong
    1676             :  * to domain */
    1677             : static GEN
    1678        4039 : lfun_OK(GEN linit, GEN s, GEN sdom, long bitprec)
    1679             : {
    1680        4039 :   GEN N, gas, S, FVga, res, ss = s;
    1681        4039 :   long prec = nbits2prec(bitprec);
    1682             : 
    1683        4039 :   FVga = lfun_get_factgammavec(linit_get_tech(linit));
    1684        4039 :   S = lfunlambda_OK(linit, s, sdom, bitprec);
    1685        4039 :   if (typ(S)==t_SER)
    1686             :   {
    1687        1365 :     long d = lg(S) - 2 + fracgammadegree(FVga);
    1688        1365 :     if (typ(s) == t_SER)
    1689         973 :       ss = sertoser(s, d);
    1690             :     else
    1691         392 :       ss = deg1ser_shallow(gen_1, s, varn(S), d);
    1692             :   }
    1693        4039 :   gas = gammafactproduct(FVga, ss, prec);
    1694        4039 :   N = ldata_get_conductor(linit_get_ldata(linit));
    1695        4039 :   res = gdiv(S, gmul(gpow(N, gdivgs(ss, 2), prec), gas));
    1696        4039 :   if (typ(s)!=t_SER && typ(res)==t_SER)
    1697             :   {
    1698         427 :     long v = valp(res);
    1699         427 :     if (v > 0) return gen_0;
    1700         392 :     if (v == 0) res = gel(res, 2);
    1701             :     else
    1702         273 :       setlg(res, minss(lg(res), 2-v));
    1703             :   }
    1704        4004 :   return gprec_w(res, prec);
    1705             : }
    1706             : 
    1707             : GEN
    1708        3248 : lfun(GEN lmisc, GEN s, long bitprec)
    1709             : {
    1710        3248 :   pari_sp av = avma;
    1711             :   GEN linit, dom, z;
    1712             :   long der;
    1713        3248 :   s = get_domain(s, &dom, &der);
    1714        3248 :   linit = lfuninit(lmisc, dom, der, bitprec);
    1715        3241 :   z = lfun_OK(linit, s, dom, bitprec);
    1716        3241 :   return gerepilecopy(av, z);
    1717             : }
    1718             : 
    1719             : /* given a t_SER a+x*s(x), return x*s(x), shallow */
    1720             : static GEN
    1721          42 : sersplit1(GEN s, GEN *head)
    1722             : {
    1723          42 :   long i, l = lg(s);
    1724             :   GEN y;
    1725          42 :   *head = simplify_shallow(mysercoeff(s, 0));
    1726          42 :   if (valp(s) > 0) return s;
    1727          28 :   y = cgetg(l-1, t_SER); y[1] = s[1];
    1728          28 :   setvalp(y, 1);
    1729          28 :   for (i=3; i < l; i++) gel(y,i-1) = gel(s,i);
    1730          28 :   return normalize(y);
    1731             : }
    1732             : 
    1733             : /* order of pole of Lambda at s (0 if regular point) */
    1734             : static long
    1735        1848 : lfunlambdaord(GEN linit, GEN s)
    1736             : {
    1737        1848 :   GEN tech = linit_get_tech(linit);
    1738        1848 :   if (linit_get_type(linit)==t_LDESC_PRODUCT)
    1739             :   {
    1740         224 :     GEN v = lfunprod_get_fact(linit_get_tech(linit));
    1741         224 :     GEN F = gel(v, 1), E = gel(v, 2), C = gel(v, 3);
    1742         224 :     long i, ex = 0, l = lg(F);
    1743         840 :     for (i = 1; i < l; i++)
    1744         616 :       ex += lfunlambdaord(gel(F,i), s) * (E[i]+C[i]);
    1745         224 :     return ex;
    1746             :   }
    1747        1624 :   if (ldata_get_residue(linit_get_ldata(linit)))
    1748             :   {
    1749         469 :     GEN r = lfunpoleresidue(lfun_get_Residue(tech), s);
    1750         469 :     if (r) return lg(r)-2;
    1751             :   }
    1752        1533 :   return 0;
    1753             : }
    1754             : 
    1755             : /* derivative of order m > 0 of L (flag = 0) or Lambda (flag = 1) */
    1756             : static GEN
    1757        1281 : lfunderiv(GEN lmisc, long m, GEN s, long flag, long bitprec)
    1758             : {
    1759        1281 :   pari_sp ltop = avma;
    1760        1281 :   GEN res, S = NULL, linit, dom;
    1761        1281 :   long der, prec = nbits2prec(bitprec);
    1762        1281 :   if (m <= 0) pari_err_DOMAIN("lfun", "D", "<=", gen_0, stoi(m));
    1763        1274 :   s = get_domain(s, &dom, &der);
    1764        1274 :   linit = lfuninit(lmisc, dom, der + m, bitprec);
    1765        1274 :   if (typ(s) == t_SER)
    1766             :   {
    1767          42 :     long v, l = lg(s)-1;
    1768             :     GEN sh;
    1769          42 :     if (valp(s) < 0) pari_err_DOMAIN("lfun","valuation", "<", gen_0, s);
    1770          42 :     S = sersplit1(s, &sh);
    1771          42 :     v = valp(S);
    1772          42 :     s = deg1ser_shallow(gen_1, sh, varn(S), m + (l+v-1)/v);
    1773             :   }
    1774             :   else
    1775             :   {
    1776        1232 :     long ex = lfunlambdaord(linit, s);
    1777             :     /* HACK: pretend lfuninit was done to right accuracy */
    1778        1232 :     s = deg1ser_shallow(gen_1, s, 0, m+1+ex);
    1779             :   }
    1780        2072 :   res = flag ? lfunlambda_OK(linit, s, dom, bitprec):
    1781         798 :                lfun_OK(linit, s, dom, bitprec);
    1782        1274 :   if (S)
    1783          42 :     res = gsubst(derivn(res, m, -1), varn(S), S);
    1784        1232 :   else if (typ(res)==t_SER)
    1785             :   {
    1786        1232 :     long v = valp(res);
    1787        1232 :     if (v > m) { set_avma(ltop); return gen_0; }
    1788        1225 :     if (v >= 0)
    1789        1155 :       res = gmul(mysercoeff(res, m), mpfact(m));
    1790             :     else
    1791          70 :       res = derivn(res, m, -1);
    1792             :   }
    1793        1267 :   return gerepilecopy(ltop, gprec_w(res, prec));
    1794             : }
    1795             : 
    1796             : GEN
    1797        1211 : lfunlambda0(GEN lmisc, GEN s, long der, long bitprec)
    1798             : {
    1799             :   return der? lfunderiv(lmisc, der, s, 1, bitprec)
    1800        1211 :             : lfunlambda(lmisc, s, bitprec);
    1801             : }
    1802             : 
    1803             : GEN
    1804        3241 : lfun0(GEN lmisc, GEN s, long der, long bitprec)
    1805             : {
    1806             :   return der? lfunderiv(lmisc, der, s, 0, bitprec)
    1807        3241 :             : lfun(lmisc, s, bitprec);
    1808             : }
    1809             : 
    1810             : GEN
    1811       11520 : lfunhardy(GEN lmisc, GEN t, long bitprec)
    1812             : {
    1813       11520 :   pari_sp ltop = avma;
    1814       11520 :   long prec = nbits2prec(bitprec), d;
    1815             :   GEN argz, z, linit, ldata, tech, dom, w2, k2, expot, h, a, k;
    1816             : 
    1817       11520 :   switch(typ(t))
    1818             :   {
    1819       11513 :     case t_INT: case t_FRAC: case t_REAL: break;
    1820           7 :     default: pari_err_TYPE("lfunhardy",t);
    1821             :   }
    1822             : 
    1823       11513 :   ldata = lfunmisc_to_ldata_shallow(lmisc);
    1824       11513 :   if (!is_linit(lmisc)) lmisc = ldata;
    1825       11513 :   k = ldata_get_k(ldata);
    1826       11513 :   d = ldata_get_degree(ldata);
    1827       11513 :   dom = mkvec3(gmul2n(k, -1), gen_0, gabs(t,LOWDEFAULTPREC));
    1828       11513 :   linit = lfuninit(lmisc, dom, 0, bitprec);
    1829       11513 :   tech = linit_get_tech(linit);
    1830       11513 :   w2 = lfun_get_w2(tech);
    1831       11513 :   k2 = lfun_get_k2(tech);
    1832       11513 :   expot = lfun_get_expot(tech);
    1833       11513 :   z = mkcomplex(k2, t);
    1834       11513 :   argz = gatan(gdiv(t, k2), prec); /* more accurate than garg since k/2 \in Q */
    1835             :   /* prec may have increased: don't lose accuracy if |z|^2 is exact */
    1836       11513 :   prec = precision(argz);
    1837       11513 :   a = gsub(gmulsg(d, gmul(t, gmul2n(argz,-1))),
    1838             :            gmul(expot,glog(gnorm(z),prec)));
    1839       11513 :   h = lfunlambda_OK(linit, z, mkvec(t), bitprec);
    1840       11513 :   if (!isint1(w2) && typ(ldata_get_dual(ldata))==t_INT)
    1841        9562 :     h = mulreal(h, w2);
    1842       11513 :   if (typ(h) == t_COMPLEX && gexpo(imag_i(h)) < -(bitprec >> 1))
    1843        1881 :     h = real_i(h);
    1844       11513 :   return gerepileupto(ltop, gmul(h, gexp(a, prec)));
    1845             : }
    1846             : 
    1847             : /* L = log(t); return  \sum_{i = 0}^{v-1}  R[-i-1] L^i/i! */
    1848             : static GEN
    1849        3787 : theta_pole_contrib(GEN R, long v, GEN L)
    1850             : {
    1851        3787 :   GEN s = mysercoeff(R,-v);
    1852             :   long i;
    1853        3948 :   for (i = v-1; i >= 1; i--)
    1854         161 :     s = gadd(mysercoeff(R,-i), gdivgs(gmul(s,L), i));
    1855        3787 :   return s;
    1856             : }
    1857             : /* subtract successively rather than adding everything then subtracting.
    1858             :  * The polar part is "large" and suffers from cancellation: a little stabler
    1859             :  * this way */
    1860             : static GEN
    1861        4333 : theta_add_polar_part(GEN S, GEN R, GEN t, long prec)
    1862             : {
    1863        4333 :   GEN logt = NULL;
    1864        4333 :   long j, l = lg(R);
    1865        8120 :   for (j = 1; j < l; j++)
    1866             :   {
    1867        3787 :     GEN Rj = gel(R,j), b = gel(Rj,1), Rb = gel(Rj,2);
    1868        3787 :     long v = -valp(Rb);
    1869        3787 :     if (v > 1 && !logt) logt = glog(t, prec);
    1870        3787 :     S = gsub(S, gmul(theta_pole_contrib(Rb,v,logt), gpow(t,b,prec)));
    1871             :   }
    1872        4333 :   return S;
    1873             : }
    1874             : 
    1875             : static long
    1876        2499 : lfuncheckfeq_i(GEN theta, GEN thetad, GEN t0, GEN t0i, long bitprec)
    1877             : {
    1878        2499 :   GEN ldata = linit_get_ldata(theta);
    1879             :   GEN S0, S0i, w, eno;
    1880        2499 :   long prec = nbits2prec(bitprec);
    1881        2499 :   if (thetad)
    1882          35 :     S0 = lfuntheta(thetad, t0, 0, bitprec);
    1883             :   else
    1884        2464 :     S0 = conj_i(lfuntheta(theta, conj_i(t0), 0, bitprec));
    1885        2499 :   S0i = lfuntheta(theta, t0i, 0, bitprec);
    1886             : 
    1887        2499 :   eno = ldata_get_rootno(ldata);
    1888        2499 :   if (ldata_get_residue(ldata))
    1889             :   {
    1890         476 :     GEN R = theta_get_R(linit_get_tech(theta));
    1891         476 :     if (gequal0(R))
    1892             :     {
    1893             :       GEN v, r;
    1894          49 :       if (ldata_get_type(ldata) == t_LFUN_NF)
    1895             :       { /* inefficient since theta not needed; no need to optimize for this
    1896             :            (artificial) query [e.g. lfuncheckfeq(t_POL)] */
    1897          21 :         GEN T = gel(ldata_get_an(ldata), 2);
    1898          21 :         GEN L = lfunzetakinit(T,zerovec(3),0,0,bitprec);
    1899          21 :         return lfuncheckfeq(L,t0,bitprec);
    1900             :       }
    1901          28 :       v = lfunrootres(theta, bitprec);
    1902          28 :       r = gel(v,1);
    1903          28 :       if (gequal0(eno)) eno = gel(v,3);
    1904          28 :       R = lfunrtoR_i(ldata, r, eno, nbits2prec(bitprec));
    1905             :     }
    1906         455 :     S0i = theta_add_polar_part(S0i, R, t0, prec);
    1907             :   }
    1908        2478 :   if (gequal0(S0i) || gequal0(S0)) pari_err_PREC("lfuncheckfeq");
    1909        2478 :   w = gdiv(S0i, gmul(S0, gpow(t0, ldata_get_k(ldata), prec)));
    1910             :   /* missing rootno: guess it */
    1911        2478 :   if (gequal0(eno)) eno = lfunrootno(theta, bitprec);
    1912        2478 :   w = gsub(w, eno);
    1913        2478 :   if (thetad) w = gdiv(w, eno); /* |eno| may be large in non-dual case */
    1914        2478 :   return gexpo(w);
    1915             : }
    1916             : 
    1917             : /* Check whether the coefficients, conductor, weight, polar part and root
    1918             :  * number are compatible with the functional equation at t0 and 1/t0.
    1919             :  * Different from lfunrootres. */
    1920             : long
    1921        2569 : lfuncheckfeq(GEN lmisc, GEN t0, long bitprec)
    1922             : {
    1923             :   GEN ldata, theta, thetad, t0i;
    1924             :   pari_sp av;
    1925             : 
    1926        2569 :   if (is_linit(lmisc) && linit_get_type(lmisc)==t_LDESC_PRODUCT)
    1927             :   {
    1928         112 :     GEN v = lfunprod_get_fact(linit_get_tech(lmisc)), F = gel(v,1);
    1929         112 :     long i, b = -bitprec, l = lg(F);
    1930         112 :     for (i = 1; i < l; i++) b = maxss(b, lfuncheckfeq(gel(F,i), t0, bitprec));
    1931         112 :     return b;
    1932             :   }
    1933        2457 :   av = avma;
    1934        2457 :   if (!t0)
    1935             :   { /* ~Pi/3 + I/7, some random complex number */
    1936        2387 :     t0 = mkcomplex(sstoQ(355,339), sstoQ(1,7));
    1937        2387 :     t0i = ginv(t0);
    1938             :   }
    1939          70 :   else if (gcmpgs(gnorm(t0), 1) < 0) { t0i = t0; t0 = ginv(t0); }
    1940          63 :   else t0i = ginv(t0);
    1941             :   /* |t0| >= 1 */
    1942        2457 :   theta = lfunthetacheckinit(lmisc, t0i, 0, bitprec);
    1943        2457 :   ldata = linit_get_ldata(theta);
    1944        2457 :   thetad = theta_dual(theta, ldata_get_dual(ldata));
    1945        2457 :   return gc_long(av, lfuncheckfeq_i(theta, thetad, t0, t0i, bitprec));
    1946             : }
    1947             : 
    1948             : /*******************************************************************/
    1949             : /*       Compute root number and residues                          */
    1950             : /*******************************************************************/
    1951             : /* round root number to \pm 1 if close to integer. */
    1952             : static GEN
    1953        3899 : ropm1(GEN eno, long prec)
    1954             : {
    1955             :   long e;
    1956        3899 :   GEN r = grndtoi(eno, &e);
    1957        3899 :   return (e < -prec2nbits(prec)/2)? r: eno;
    1958             : }
    1959             : 
    1960             : /* theta for t=1/sqrt(2) and t2==2t simultaneously, saving 25% of the work.
    1961             :  * Assume correct initialization (no thetacheck) */
    1962             : static void
    1963          91 : lfunthetaspec(GEN linit, long bitprec, GEN *pv, GEN *pv2)
    1964             : {
    1965          91 :   pari_sp av = avma, av2;
    1966             :   GEN t, Vga, an, K, ldata, thetainit, v, v2, vroots;
    1967             :   long L, prec, n, d;
    1968             : 
    1969          91 :   ldata = linit_get_ldata(linit);
    1970          91 :   thetainit = linit_get_tech(linit);
    1971          91 :   prec = nbits2prec(bitprec);
    1972          91 :   Vga = ldata_get_gammavec(ldata); d = lg(Vga)-1;
    1973          91 :   if (vgaeasytheta(Vga))
    1974             :   {
    1975          70 :     GEN v2 = sqrtr(real2n(1, nbits2prec(bitprec)));
    1976          70 :     GEN v = shiftr(v2,-1);
    1977          70 :     *pv = lfuntheta(linit, v,  0, bitprec);
    1978          70 :     *pv2= lfuntheta(linit, v2, 0, bitprec);
    1979          70 :     return;
    1980             :   }
    1981          21 :   an = RgV_kill0( theta_get_an(thetainit) );
    1982          21 :   L = lg(an)-1;
    1983             :   /* to compute theta(1/sqrt(2)) */
    1984          21 :   t = ginv(gsqrt(gmul2n(ldata_get_conductor(ldata), 1), prec));
    1985             :   /* t = 1/sqrt(2N) */
    1986             : 
    1987             :   /* From then on, the code is generic and could be used to compute
    1988             :    * theta(t) / theta(2t) without assuming t = 1/sqrt(2) */
    1989          21 :   K = theta_get_K(thetainit);
    1990          21 :   vroots = mkvroots(d, L, prec);
    1991          21 :   t = gpow(t, gdivgs(gen_2, d), prec); /* rt variant: t->t^(2/d) */
    1992             :   /* v = \sum_{n <= L, n odd} a_n K(nt) */
    1993       83342 :   for (v = gen_0, n = 1; n <= L; n+=2)
    1994             :   {
    1995       83321 :     GEN tn, Kn, a = gel(an, n);
    1996             : 
    1997       83321 :     if (!a) continue;
    1998       12299 :     av2 = avma;
    1999       12299 :     tn = gmul(t, gel(vroots,n));
    2000       12299 :     Kn = gammamellininvrt(K, tn, bitprec);
    2001       12299 :     v = gerepileupto(av2, gadd(v, gmul(a,Kn)));
    2002             :   }
    2003             :   /* v += \sum_{n <= L, n even} a_n K(nt), v2 = \sum_{n <= L/2} a_n K(2n t) */
    2004       83328 :   for (v2 = gen_0, n = 1; n <= L/2; n++)
    2005             :   {
    2006       83307 :     GEN t2n, K2n, a = gel(an, n), a2 = gel(an,2*n);
    2007             : 
    2008       83307 :     if (!a && !a2) continue;
    2009        9198 :     av2 = avma;
    2010        9198 :     t2n = gmul(t, gel(vroots,2*n));
    2011        9198 :     K2n = gerepileupto(av2, gammamellininvrt(K, t2n, bitprec));
    2012        9198 :     if (a) v2 = gadd(v2, gmul(a, K2n));
    2013        9198 :     if (a2) v = gadd(v,  gmul(a2,K2n));
    2014             :   }
    2015          21 :   *pv = v;
    2016          21 :   *pv2 = v2;
    2017          21 :   gerepileall(av, 2, pv,pv2);
    2018             : }
    2019             : 
    2020             : static GEN
    2021          56 : Rtor(GEN a, GEN R, GEN ldata, long prec)
    2022             : {
    2023          56 :   GEN FVga = gammafactor(ldata_get_gammavec(ldata));
    2024          56 :   GEN Na = gpow(ldata_get_conductor(ldata), gdivgs(a,2), prec);
    2025          56 :   return gdiv(R, gmul(Na, gammafactproduct(FVga, a, prec)));
    2026             : }
    2027             : 
    2028             : /* v = theta~(t), vi = theta(1/t) */
    2029             : static GEN
    2030        3878 : get_eno(GEN R, GEN k, GEN t, GEN v, GEN vi, long vx, long bitprec, long force)
    2031             : {
    2032        3878 :   long prec = nbits2prec(bitprec);
    2033        3878 :   GEN a0, a1, S = deg1pol(gmul(gpow(t,k,prec), gneg(v)), vi, vx);
    2034             : 
    2035        3878 :   S = theta_add_polar_part(S, R, t, prec);
    2036        3878 :   if (typ(S) != t_POL || degpol(S) != 1) return NULL;
    2037        3878 :   a1 = gel(S,3); if (!force && gexpo(a1) < -bitprec/4) return NULL;
    2038        3843 :   a0 = gel(S,2);
    2039        3843 :   return gdiv(a0, gneg(a1));
    2040             : 
    2041             : }
    2042             : /* Return w using theta(1/t) - w t^k \bar{theta}(t) = polar_part(t,w).
    2043             :  * The full Taylor development of L must be known */
    2044             : GEN
    2045        3843 : lfunrootno(GEN linit, long bitprec)
    2046             : {
    2047             :   GEN ldata, t, eno, v, vi, R, thetad;
    2048        3843 :   long c = 0, prec = nbits2prec(bitprec), vx = fetch_var();
    2049             :   GEN k;
    2050             :   pari_sp av;
    2051             : 
    2052             :   /* initialize for t > 1/sqrt(2) */
    2053        3843 :   linit = lfunthetacheckinit(linit, dbltor(sqrt(0.5)), 0, bitprec);
    2054        3843 :   ldata = linit_get_ldata(linit);
    2055        3843 :   k = ldata_get_k(ldata);
    2056        9065 :   R = ldata_get_residue(ldata)? lfunrtoR_eno(ldata, pol_x(vx), prec)
    2057        5222 :                               : cgetg(1, t_VEC);
    2058        3843 :   t = gen_1;
    2059        3843 :   v = lfuntheta(linit, t, 0, bitprec);
    2060        3843 :   thetad = theta_dual(linit, ldata_get_dual(ldata));
    2061        3843 :   vi = !thetad ? conj_i(v): lfuntheta(thetad, t, 0, bitprec);
    2062        3843 :   eno = get_eno(R,k,t,vi,v, vx, bitprec, 0);
    2063        3843 :   if (!eno && !thetad)
    2064             :   { /* t = sqrt(2), vi = theta(1/t), v = theta(t) */
    2065          35 :     lfunthetaspec(linit, bitprec, &vi, &v);
    2066          35 :     t = sqrtr(utor(2, prec));
    2067          35 :     eno = get_eno(R,k,t,conj_i(v),vi, vx, bitprec, 0);
    2068             :   }
    2069        3843 :   av = avma;
    2070        7686 :   while (!eno)
    2071             :   {
    2072           0 :     t = addsr(1, shiftr(utor(pari_rand(), prec), -2-BITS_IN_LONG));
    2073             :     /* t in [1,1.25[ */
    2074           0 :     v = thetad? lfuntheta(thetad, t, 0, bitprec)
    2075           0 :               : conj_i(lfuntheta(linit, t, 0, bitprec));
    2076           0 :     vi = lfuntheta(linit, ginv(t), 0, bitprec);
    2077           0 :     eno = get_eno(R,k,t,v,vi, vx, bitprec, c++ == 5);
    2078           0 :     set_avma(av);
    2079             :   }
    2080        3843 :   delete_var(); return ropm1(eno,prec);
    2081             : }
    2082             : 
    2083             : /* Find root number and/or residues when L-function coefficients and
    2084             :    conductor are known. For the moment at most a single residue allowed. */
    2085             : GEN
    2086        2527 : lfunrootres(GEN data, long bitprec)
    2087             : {
    2088        2527 :   pari_sp ltop = avma;
    2089             :   GEN k, w, r, R, a, b, e, v, v2, be, ldata, linit;
    2090             :   long prec;
    2091             : 
    2092        2527 :   ldata = lfunmisc_to_ldata_shallow(data);
    2093        2527 :   r = ldata_get_residue(ldata);
    2094        2527 :   k = ldata_get_k(ldata);
    2095        2527 :   if (r) r = normalize_simple_pole(r, k);
    2096        2527 :   if (!r || residues_known(r))
    2097             :   {
    2098        2471 :     w = lfunrootno(data, bitprec);
    2099        2471 :     if (!r)
    2100        1092 :       r = R = gen_0;
    2101             :     else
    2102        1379 :       R = lfunrtoR_eno(ldata, w, nbits2prec(bitprec));
    2103        2471 :     return gerepilecopy(ltop, mkvec3(r, R, w));
    2104             :   }
    2105          56 :   linit = lfunthetacheckinit(data, dbltor(sqrt(0.5)), 0, bitprec);
    2106          56 :   prec = nbits2prec(bitprec);
    2107          56 :   if (lg(r) > 2) pari_err_IMPL("multiple poles in lfunrootres");
    2108             :   /* Now residue unknown, and r = [[be,0]]. */
    2109          56 :   be = gmael(r, 1, 1);
    2110          56 :   w = ldata_get_rootno(ldata);
    2111          56 :   if (ldata_isreal(ldata) && gequalm1(w))
    2112           0 :     R = lfuntheta(linit, gen_1, 0, bitprec);
    2113             :   else
    2114             :   {
    2115          56 :     lfunthetaspec(linit, bitprec, &v2, &v);
    2116          56 :     if (gequal(gmulsg(2, be), k)) pari_err_IMPL("pole at k/2 in lfunrootres");
    2117          56 :     if (gequal(be, k))
    2118             :     {
    2119           7 :       GEN p2k = gpow(gen_2,k,prec);
    2120           7 :       a = conj_i(gsub(gmul(p2k, v), v2));
    2121           7 :       b = subiu(p2k, 1);
    2122           7 :       e = gmul(gsqrt(p2k, prec), gsub(v2, v));
    2123             :     }
    2124             :     else
    2125             :     {
    2126          49 :       GEN p2k = gpow(gen_2,k,prec);
    2127          49 :       GEN tk2 = gsqrt(p2k, prec);
    2128          49 :       GEN tbe = gpow(gen_2, be, prec);
    2129          49 :       GEN tkbe = gpow(gen_2, gdivgs(gsub(k, be), 2), prec);
    2130          49 :       a = conj_i(gsub(gmul(tbe, v), v2));
    2131          49 :       b = gsub(gdiv(tbe, tkbe), tkbe);
    2132          49 :       e = gsub(gmul(gdiv(tbe, tk2), v2), gmul(tk2, v));
    2133             :     }
    2134          56 :     if (!isintzero(w)) R = gdiv(gsub(e, gmul(a, w)), b);
    2135             :     else
    2136             :     { /* Now residue unknown, r = [[be,0]], and w unknown. */
    2137           0 :       GEN t0  = mkfrac(stoi(11),stoi(10));
    2138           0 :       GEN th1 = lfuntheta(linit, t0,  0, bitprec);
    2139           0 :       GEN th2 = lfuntheta(linit, ginv(t0), 0, bitprec);
    2140           0 :       GEN tbe = gpow(t0, gmulsg(2, be), prec);
    2141           0 :       GEN tkbe = gpow(t0, gsub(k, be), prec);
    2142           0 :       GEN tk2 = gpow(t0, k, prec);
    2143           0 :       GEN c = conj_i(gsub(gmul(tbe, th1), th2));
    2144           0 :       GEN d = gsub(gdiv(tbe, tkbe), tkbe);
    2145           0 :       GEN f = gsub(gmul(gdiv(tbe, tk2), th2), gmul(tk2, th1));
    2146           0 :       GEN D = gsub(gmul(a, d), gmul(b, c));
    2147           0 :       w = gdiv(gsub(gmul(d, e), gmul(b, f)), D);
    2148           0 :       R = gdiv(gsub(gmul(a, f), gmul(c, e)), D);
    2149             :     }
    2150             :   }
    2151          56 :   r = normalize_simple_pole(Rtor(be, R, ldata, prec), be);
    2152          56 :   R = lfunrtoR_i(ldata, r, w, prec);
    2153          56 :   return gerepilecopy(ltop, mkvec3(r, R, ropm1(w, prec)));
    2154             : }
    2155             : 
    2156             : /*******************************************************************/
    2157             : /*                           Zeros                                 */
    2158             : /*******************************************************************/
    2159             : struct lhardyz_t {
    2160             :   long bitprec, prec;
    2161             :   GEN linit;
    2162             : };
    2163             : 
    2164             : static GEN
    2165       11023 : lfunhardyzeros(void *E, GEN t)
    2166             : {
    2167       11023 :   struct lhardyz_t *S = (struct lhardyz_t*)E;
    2168       11023 :   long prec = S->prec;
    2169       11023 :   GEN h = lfunhardy(S->linit, t, S->bitprec);
    2170       11023 :   if (typ(h) == t_REAL && realprec(h) < prec) h = gprec_w(h, prec);
    2171       11023 :   return h;
    2172             : }
    2173             : 
    2174             : /* initialize for computation on critical line up to height h, zero
    2175             :  * of order <= m */
    2176             : static GEN
    2177         420 : lfuncenterinit(GEN lmisc, double h, long m, long bitprec)
    2178             : {
    2179         420 :   if (m < 0)
    2180             :   { /* choose a sensible default */
    2181         420 :     if (!is_linit(lmisc) || linit_get_type(lmisc) != t_LDESC_INIT) m = 4;
    2182             :     else
    2183             :     {
    2184         378 :       GEN domain = lfun_get_domain(linit_get_tech(lmisc));
    2185         378 :       m = domain_get_der(domain);
    2186             :     }
    2187             :   }
    2188         420 :   return lfuninit(lmisc, mkvec(dbltor(h)), m, bitprec);
    2189             : }
    2190             : 
    2191             : long
    2192         434 : lfunorderzero(GEN lmisc, long m, long bitprec)
    2193             : {
    2194         434 :   pari_sp ltop = avma;
    2195             :   GEN eno, ldata, linit, k2;
    2196             :   long G, c0, c, st;
    2197             : 
    2198         434 :   if (is_linit(lmisc) && linit_get_type(lmisc) == t_LDESC_PRODUCT)
    2199             :   {
    2200          63 :     GEN M = gmael(linit_get_tech(lmisc), 2,1);
    2201             :     long i;
    2202          63 :     for (c=0,i=1; i < lg(M); i++) c += lfunorderzero(gel(M,i), m, bitprec);
    2203          63 :     return c;
    2204             :   }
    2205         371 :   linit = lfuncenterinit(lmisc, 0, m, bitprec);
    2206         371 :   ldata = linit_get_ldata(linit);
    2207         371 :   eno = ldata_get_rootno(ldata);
    2208         371 :   G = -bitprec/2;
    2209         371 :   c0 = 0; st = 1;
    2210         371 :   if (ldata_isreal(ldata))
    2211             :   {
    2212         308 :     if (!gequal1(eno)) c0 = 1;
    2213         308 :     st = 2;
    2214             :   }
    2215         371 :   k2 = gmul2n(ldata_get_k(ldata), -1);
    2216         392 :   for (c = c0;; c += st)
    2217         413 :     if (gexpo(lfun0(linit, k2, c, bitprec)) > G) return gc_long(ltop, c);
    2218             : }
    2219             : 
    2220             : GEN
    2221          49 : lfunzeros(GEN ldata, GEN lim, long divz, long bitprec)
    2222             : {
    2223          49 :   pari_sp ltop = avma;
    2224             :   GEN ldataf, linit, N, pi2, cN, pi2div, w, T, Vga, h1, h2;
    2225          49 :   long i, d, W, NEWD, precinit, ct, s, prec = nbits2prec(bitprec);
    2226             :   double maxt;
    2227             :   GEN maxtr, maxtr1;
    2228             :   struct lhardyz_t S;
    2229             : 
    2230          49 :   if (typ(lim) == t_VEC)
    2231             :   {
    2232          14 :     if (lg(lim) != 3 || gcmp(gel(lim,1),gel(lim,2)) >= 0
    2233          14 :                      || gcmp(gel(lim,1),gen_0) < 0)
    2234           0 :       pari_err_TYPE("lfunzeros",lim);
    2235          14 :     h1 = gel(lim,1); h2 = gel(lim,2);
    2236             :   }
    2237             :   else
    2238             :   {
    2239          35 :     if (gcmp(lim,gen_0) <= 0)
    2240           0 :       pari_err_TYPE("lfunzeros",lim);
    2241          35 :     h1 = gen_0; h2 = lim;
    2242             :   }
    2243          49 :   maxt = gtodouble(h2);
    2244             : 
    2245          49 :   if (is_linit(ldata) && linit_get_type(ldata) == t_LDESC_PRODUCT)
    2246             :   {
    2247           0 :     GEN v, M = gmael(linit_get_tech(ldata), 2,1);
    2248           0 :     long l = lg(M);
    2249           0 :     v = cgetg(l, t_VEC);
    2250           0 :     for (i = 1; i < l; i++)
    2251           0 :       gel(v,i) = lfunzeros(gel(M,i), lim, divz, bitprec);
    2252           0 :     return gerepileupto(ltop, vecsort0(shallowconcat1(v), NULL, 0));
    2253             :   }
    2254          49 :   S.linit = linit = lfuncenterinit(ldata, maxt + 1, -1, bitprec);
    2255          49 :   S.bitprec = bitprec;
    2256          49 :   S.prec = prec;
    2257          49 :   ldataf = linit_get_ldata(linit);
    2258          49 :   Vga = ldata_get_gammavec(ldataf); d = lg(Vga) - 1;
    2259          49 :   N = ldata_get_conductor(ldataf);
    2260          49 :   NEWD = minss((long) ceil(bitprec+(M_PI/(4*M_LN2))*d*maxt),
    2261             :                lfun_get_bitprec(linit_get_tech(linit)));
    2262          49 :   precinit = prec; prec = nbits2prec(NEWD);
    2263          49 :   pi2 = Pi2n(1, prec);
    2264          49 :   cN = gdiv(N, gpowgs(Pi2n(-1, prec), d));
    2265          49 :   cN = gexpo(cN) >= 0? gaddsg(d, gmulsg(2, glog(cN, prec))): stoi(d);
    2266          49 :   pi2div = gdivgs(pi2, labs(divz));
    2267          49 :   ct = 0;
    2268          49 :   T = h1;
    2269          49 :   if (gequal0(h1))
    2270             :   {
    2271          42 :     GEN r = ldata_get_residue(ldataf);
    2272          42 :     if (!r || gequal0(r))
    2273             :     {
    2274          28 :       ct = lfunorderzero(linit, -1, bitprec);
    2275          28 :       if (ct) T = real2n(-prec2nbits(prec)/(2*ct), prec);
    2276             :     }
    2277             :   }
    2278             :   /* initialize for 100 further zeros, double later if needed */
    2279          49 :   W = 100 + ct; w = cgetg(W+1,t_VEC);
    2280          49 :   for (i=1; i<=ct; i++) gel(w,i) = gen_0;
    2281          49 :   s = gsigne(lfunhardyzeros(&S, T));
    2282          49 :   maxtr = h2; maxtr1 = gaddsg(1, maxtr);
    2283         448 :   while (gcmp(T, maxtr1) < 0)
    2284             :   {
    2285         399 :     pari_sp av = avma;
    2286         399 :     GEN T0 = T, z;
    2287             :     for(;;)
    2288        6195 :     {
    2289             :       long s0;
    2290             :       GEN L;
    2291        6594 :       if (gcmp(T, pi2) >= 0)
    2292        4221 :         L = gadd(cN, gmulsg(d, glog(gdiv(T, pi2), prec)));
    2293             :       else
    2294        2373 :         L = cN;
    2295        6594 :       T = gadd(T, gdiv(pi2div, L));
    2296        6594 :       if (gcmp(T, maxtr1) > 0) goto END;
    2297        6573 :       s0 = gsigne(lfunhardyzeros(&S, T));
    2298        6573 :       if (s0 != s) { s = s0; break; }
    2299             :     }
    2300         378 :     T = gerepileupto(av, T);
    2301         378 :     z = zbrent(&S, lfunhardyzeros, T0, T, prec);
    2302         378 :     if (gcmp(z, maxtr) > 0) break;
    2303         350 :     if (typ(z) == t_REAL) z  = rtor(z, precinit);
    2304             :     /* room for twice as many zeros */
    2305         350 :     if (ct >= W) { W *= 2; w = vec_lengthen(w, W); }
    2306         350 :     gel(w, ++ct) = z;
    2307             :   }
    2308             : END:
    2309          49 :   setlg(w, ct+1); return gerepilecopy(ltop, w);
    2310             : }
    2311             : 
    2312             : /*******************************************************************/
    2313             : /*       Guess conductor                                           */
    2314             : /*******************************************************************/
    2315             : struct huntcond_t {
    2316             :   GEN k;
    2317             :   GEN theta, thetad;
    2318             :   GEN *pM, *psqrtM, *pMd, *psqrtMd;
    2319             : };
    2320             : 
    2321             : static void
    2322       10789 : condset(struct huntcond_t *S, GEN M, long prec)
    2323             : {
    2324       10789 :   *(S->pM) = M;
    2325       10789 :   *(S->psqrtM) = gsqrt(ginv(M), prec);
    2326       10789 :   if (S->thetad != S->theta)
    2327             :   {
    2328           0 :     *(S->pMd) = *(S->pM);
    2329           0 :     *(S->psqrtMd) = *(S->psqrtM);
    2330             :   }
    2331       10789 : }
    2332             : 
    2333             : /* M should eventually converge to N, the conductor. L has no pole. */
    2334             : static GEN
    2335        6468 : wrap1(void *E, GEN M)
    2336             : {
    2337        6468 :   struct huntcond_t *S = (struct huntcond_t*)E;
    2338             :   GEN thetainit, tk, p1, p1inv;
    2339        6468 :   GEN t = mkfrac(stoi(11), stoi(10));
    2340             :   long prec, bitprec;
    2341             : 
    2342        6468 :   thetainit = linit_get_tech(S->theta);
    2343        6468 :   bitprec = theta_get_bitprec(thetainit);
    2344        6468 :   prec = nbits2prec(bitprec);
    2345        6468 :   condset(S, M, prec);
    2346        6468 :   tk = gpow(t, S->k, prec);
    2347        6468 :   p1 = lfuntheta(S->thetad, t, 0, bitprec);
    2348        6468 :   p1inv = lfuntheta(S->theta, ginv(t), 0, bitprec);
    2349        6468 :   return glog(gabs(gmul(tk, gdiv(p1, p1inv)), prec), prec);
    2350             : }
    2351             : 
    2352             : /* M should eventually converge to N, the conductor. L has a pole. */
    2353             : static GEN
    2354        4279 : wrap2(void *E, GEN M)
    2355             : {
    2356        4279 :   struct huntcond_t *S = (struct huntcond_t*)E;
    2357             :   GEN t1k, t2k, p1, p1inv, p2, p2inv, thetainit, R;
    2358        4279 :   GEN t1 = mkfrac(stoi(11), stoi(10)), t2 = mkfrac(stoi(13), stoi(11));
    2359             :   GEN t1be, t2be, t1bemk, t2bemk, t1kmbe, t2kmbe;
    2360             :   GEN F11, F12, F21, F22, P1, P2, res;
    2361             :   long prec, bitprec;
    2362        4279 :   GEN k = S->k;
    2363             : 
    2364        4279 :   thetainit = linit_get_tech(S->theta);
    2365        4279 :   bitprec = theta_get_bitprec(thetainit);
    2366        4279 :   prec = nbits2prec(bitprec);
    2367        4279 :   condset(S, M, prec);
    2368             : 
    2369        4279 :   p1 = lfuntheta(S->thetad, t1, 0, bitprec);
    2370        4279 :   p2 = lfuntheta(S->thetad, t2, 0, bitprec);
    2371        4279 :   p1inv = lfuntheta(S->theta, ginv(t1), 0, bitprec);
    2372        4279 :   p2inv = lfuntheta(S->theta, ginv(t2), 0, bitprec);
    2373        4279 :   t1k = gpow(t1, k, prec);
    2374        4279 :   t2k = gpow(t2, k, prec);
    2375        4279 :   R = theta_get_R(thetainit);
    2376        4279 :   if (typ(R) == t_VEC)
    2377             :   {
    2378        4279 :     GEN be = gmael(R, 1, 1);
    2379        4279 :     t1be = gpow(t1, be, prec); t1bemk = gdiv(gsqr(t1be), t1k);
    2380        4279 :     t2be = gpow(t2, be, prec); t2bemk = gdiv(gsqr(t2be), t2k);
    2381        4279 :     t1kmbe = gdiv(t1k, t1be);
    2382        4279 :     t2kmbe = gdiv(t2k, t2be);
    2383             :   }
    2384             :   else
    2385             :   { /* be = k */
    2386           0 :     t1be = t1k; t1bemk = t1k; t1kmbe = gen_1;
    2387           0 :     t2be = t2k; t2bemk = t2k; t2kmbe = gen_1;
    2388             :   }
    2389        4279 :   F11 = conj_i(gsub(gmul(gsqr(t1be), p1), p1inv));
    2390        4279 :   F12 = conj_i(gsub(gmul(gsqr(t2be), p2), p2inv));
    2391        4279 :   F21 = gsub(gmul(t1k, p1), gmul(t1bemk, p1inv));
    2392        4279 :   F22 = gsub(gmul(t2k, p2), gmul(t2bemk, p2inv));
    2393        4279 :   P1 = gsub(gmul(t1bemk, t1be), t1kmbe);
    2394        4279 :   P2 = gsub(gmul(t2bemk, t2be), t2kmbe);
    2395        4279 :   res = gdiv(gsub(gmul(P2,F21), gmul(P1,F22)),
    2396             :              gsub(gmul(P2,F11), gmul(P1,F12)));
    2397        4279 :   return glog(gabs(res, prec), prec);
    2398             : }
    2399             : 
    2400             : /* If flag = 0 (default) return all conductors found as integers. If
    2401             : flag = 1, return the approximations, not the integers. If flag = 2,
    2402             : return all, even nonintegers. */
    2403             : 
    2404             : static GEN
    2405          84 : checkconductor(GEN v, long bit, long flag)
    2406             : {
    2407             :   GEN w;
    2408          84 :   long e, j, k, l = lg(v);
    2409          84 :   if (flag == 2) return v;
    2410          84 :   w = cgetg(l, t_VEC);
    2411         322 :   for (j = k = 1; j < l; j++)
    2412             :   {
    2413         238 :     GEN N = grndtoi(gel(v,j), &e);
    2414         238 :     if (e < -bit) gel(w,k++) = flag ? gel(v,j): N;
    2415             :   }
    2416          84 :   if (k == 2) return gel(w,1);
    2417           7 :   setlg(w,k); return w;
    2418             : }
    2419             : 
    2420             : static GEN
    2421          91 : parse_maxcond(GEN maxN)
    2422             : {
    2423             :   GEN M;
    2424          91 :   if (!maxN)
    2425          49 :     M = utoipos(10000);
    2426          42 :   else if (typ(maxN) == t_VEC)
    2427             :   {
    2428           7 :     if (!RgV_is_ZV(maxN)) pari_err_TYPE("lfunconductor",maxN);
    2429           7 :     return ZV_sort(maxN);
    2430             :   }
    2431             :   else
    2432          35 :     M = maxN;
    2433          84 :   return (typ(M) == t_INT)? addiu(M, 1): gceil(M);
    2434             : }
    2435             : 
    2436             : GEN
    2437          91 : lfunconductor(GEN data, GEN maxcond, long flag, long bitprec)
    2438             : {
    2439             :   struct huntcond_t S;
    2440          91 :   pari_sp av = avma;
    2441          91 :   GEN ldata = lfunmisc_to_ldata_shallow(data);
    2442          91 :   GEN ld, r, v, theta, thetad, M, tdom, t0 = NULL, t0i = NULL;
    2443             :   GEN (*eval)(void *, GEN);
    2444             :   long prec;
    2445          91 :   M = parse_maxcond(maxcond);
    2446          91 :   r = ldata_get_residue(ldata);
    2447          91 :   if (typ(M) == t_VEC) /* select in list */
    2448           7 :   { eval = NULL; tdom = dbltor(0.7); }
    2449          84 :   else if (!r) { eval = wrap1; tdom = sstoQ(10,11); }
    2450             :   else
    2451             :   {
    2452          21 :     if (typ(r) == t_VEC && lg(r) > 2)
    2453           0 :       pari_err_IMPL("multiple poles in lfunconductor");
    2454          21 :     eval = wrap2; tdom = sstoQ(11,13);
    2455             :   }
    2456          91 :   if (eval) bitprec += bitprec/2;
    2457          91 :   prec = nbits2prec(bitprec);
    2458          91 :   ld = shallowcopy(ldata);
    2459          91 :   gel(ld, 5) = eval? M: gel(M,lg(M)-1);
    2460          91 :   theta = lfunthetainit_i(ld, tdom, 0, bitprec);
    2461          91 :   thetad = theta_dual(theta, ldata_get_dual(ldata));
    2462          91 :   gel(theta,3) = shallowcopy(linit_get_tech(theta));
    2463          91 :   S.k = ldata_get_k(ldata);
    2464          91 :   S.theta = theta;
    2465          91 :   S.thetad = thetad? thetad: theta;
    2466          91 :   S.pM = &gel(linit_get_ldata(theta),5);
    2467          91 :   S.psqrtM = &gel(linit_get_tech(theta),7);
    2468          91 :   if (thetad)
    2469             :   {
    2470           0 :     S.pMd = &gel(linit_get_ldata(thetad),5);
    2471           0 :     S.psqrtMd = &gel(linit_get_tech(thetad),7);
    2472             :   }
    2473          91 :   if (!eval)
    2474             :   {
    2475           7 :     long i, besti = 0, beste = -10, l = lg(M);
    2476           7 :     t0 = sstoQ(11,10); t0i = sstoQ(10,11);
    2477          49 :     for (i = 1; i < l; i++)
    2478             :     {
    2479          42 :       pari_sp av2 = avma;
    2480             :       long e;
    2481          42 :       condset(&S, gel(M,i), prec);
    2482          42 :       e = lfuncheckfeq_i(theta, thetad, t0, t0i, bitprec);
    2483          42 :       set_avma(av2);
    2484          42 :       if (e < beste) { beste = e; besti = i; }
    2485          35 :       else if (e == beste) beste = besti = 0; /* tie: forget */
    2486             :     }
    2487           7 :     if (!besti) { set_avma(av); return cgetg(1,t_VEC); }
    2488           7 :     return gerepilecopy(av, mkvec2(gel(M,besti), stoi(beste)));
    2489             :   }
    2490          84 :   v = solvestep((void*)&S, eval, ghalf, M, gen_2, 14, prec);
    2491          84 :   return gerepilecopy(av, checkconductor(v, bitprec/2, flag));
    2492             : }
    2493             : 
    2494             : /* assume chi primitive */
    2495             : static GEN
    2496         602 : znchargauss_i(GEN G, GEN chi, long bitprec)
    2497             : {
    2498         602 :   GEN z, q, F = znstar_get_N(G);
    2499             :   long prec;
    2500             : 
    2501         602 :   if (equali1(F)) return gen_1;
    2502         357 :   prec = nbits2prec(bitprec);
    2503         357 :   q = sqrtr_abs(itor(F, prec));
    2504         357 :   z = lfuntheta(mkvec2(G,chi), gen_1, 0, bitprec);
    2505         357 :   if (gexpo(z) < 10 - bitprec)
    2506             :   {
    2507          28 :     if (equaliu(F,300))
    2508             :     {
    2509          14 :       GEN z = rootsof1u_cx(25, prec);
    2510          14 :       GEN n = znconreyexp(G, chi);
    2511          14 :       if (equaliu(n, 131)) return gmul(q, gpowgs(z,14));
    2512           7 :       if (equaliu(n, 71)) return gmul(q, gpowgs(z,11));
    2513             :     }
    2514          14 :     if (equaliu(F,600))
    2515             :     {
    2516          14 :       GEN z = rootsof1u_cx(25, prec);
    2517          14 :       GEN n = znconreyexp(G, chi);
    2518          14 :       if (equaliu(n, 491)) return gmul(q, gpowgs(z,7));
    2519           7 :       if (equaliu(n, 11)) return gmul(q, gpowgs(z,18));
    2520             :     }
    2521           0 :     pari_err_BUG("znchargauss [ Theta(chi,1) = 0 ]");
    2522             :   }
    2523         329 :   z = gmul(gdiv(z, conj_i(z)), q);
    2524         329 :   if (zncharisodd(G,chi)) z = mulcxI(z);
    2525         329 :   return z;
    2526             : }
    2527             : static GEN
    2528         602 : Z_radical(GEN N, long *om)
    2529             : {
    2530         602 :   GEN P = gel(Z_factor(N), 1);
    2531         602 :   *om = lg(P)-1; return ZV_prod(P);
    2532             : }
    2533             : GEN
    2534         819 : znchargauss(GEN G, GEN chi, GEN a, long bitprec)
    2535             : {
    2536             :   GEN v, T, N, F, b0, b1, b2, bF, a1, aF, A, r, GF, tau, B, faB, u, S;
    2537         819 :   long omb0, prec = nbits2prec(bitprec);
    2538         819 :   pari_sp av = avma;
    2539             : 
    2540         819 :   if (typ(chi) != t_COL) chi = znconreylog(G,chi);
    2541         819 :   T = znchartoprimitive(G, chi);
    2542         819 :   GF  = gel(T,1);
    2543         819 :   chi = gel(T,2); /* now primitive */
    2544         819 :   N = znstar_get_N(G);
    2545         819 :   F = znstar_get_N(GF);
    2546         819 :   if (equalii(N,F)) b1 = bF = gen_1;
    2547             :   else
    2548             :   {
    2549         231 :     v = Z_ppio(diviiexact(N,F), F);
    2550         231 :     bF = gel(v,2); /* (N/F, F^oo) */
    2551         231 :     b1 = gel(v,3); /* cofactor */
    2552             :   }
    2553         819 :   if (!a) a = a1 = aF = gen_1;
    2554             :   else
    2555             :   {
    2556         770 :     if (typ(a) != t_INT) pari_err_TYPE("znchargauss",a);
    2557         770 :     a = modii(a, N);
    2558         770 :     v = Z_ppio(a, F);
    2559         770 :     aF = gel(v,2);
    2560         770 :     a1 = gel(v,3);
    2561             :   }
    2562         819 :   if (!equalii(aF, bF)) { set_avma(av); return gen_0; }
    2563         602 :   b0 = Z_radical(b1, &omb0);
    2564         602 :   b2 = diviiexact(b1, b0);
    2565         602 :   A = dvmdii(a1, b2, &r);
    2566         602 :   if (r != gen_0) { set_avma(av); return gen_0; }
    2567         602 :   B = gcdii(A,b0); faB = Z_factor(B); /* squarefree */
    2568         602 :   S = eulerphi(mkvec2(B,faB));
    2569         602 :   if (odd(omb0 + lg(gel(faB,1))-1)) S = negi(S); /* moebius(b0/B) * phi(B) */
    2570         602 :   S = mulii(S, mulii(aF,b2));
    2571         602 :   tau = znchargauss_i(GF, chi, bitprec);
    2572         602 :   u = Fp_div(b0, A, F);
    2573         602 :   if (!equali1(u))
    2574             :   {
    2575         252 :     GEN ord = zncharorder(GF, chi), z = rootsof1_cx(ord, prec);
    2576         252 :     tau = gmul(tau, znchareval(GF, chi, u, mkvec2(z,ord)));
    2577             :   }
    2578         602 :   return gerepileupto(av, gmul(tau, S));
    2579             : }

Generated by: LCOV version 1.13