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

Generated by: LCOV version 1.13