Code coverage tests

This page documents the degree to which the PARI/GP source code is tested by our public test suite, distributed with the source distribution in directory src/test/. This is measured by the gcov utility; we then process gcov output using the lcov frond-end.

We test a few variants depending on Configure flags on the pari.math.u-bordeaux.fr machine (x86_64 architecture), and agregate them in the final report:

The target is to exceed 90% coverage for all mathematical modules (given that branches depending on DEBUGLEVEL or DEBUGMEM are not covered). This script is run to produce the results below.

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

Generated by: LCOV version 1.13