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

Generated by: LCOV version 1.11