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.18.1 lcov report (development 30489-4d12223a6e) Lines: 1595 1636 97.5 %
Date: 2025-09-14 09:22:35 Functions: 163 163 100.0 %
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; either version 2 of the License, or (at your option) any later
       8             : version. It is distributed in the hope that it will be useful, but WITHOUT
       9             : ANY WARRANTY WHATSOEVER.
      10             : 
      11             : Check the License for details. You should have received a copy of it, along
      12             : with the package; see the file 'COPYING'. If not, write to the Free Software
      13             : Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. */
      14             : 
      15             : /********************************************************************/
      16             : /**                                                                **/
      17             : /**                       L-functions                              **/
      18             : /**                                                                **/
      19             : /********************************************************************/
      20             : 
      21             : #include "pari.h"
      22             : #include "paripriv.h"
      23             : 
      24             : #define DEBUGLEVEL DEBUGLEVEL_lfun
      25             : 
      26             : /*******************************************************************/
      27             : /*  Accessors                                                      */
      28             : /*******************************************************************/
      29             : 
      30             : static GEN
      31       12017 : mysercoeff(GEN x, long n)
      32             : {
      33       12017 :   long N = n - valser(x);
      34       12017 :   return (N < 0)? gen_0: gel(x, N+2);
      35             : }
      36             : 
      37             : long
      38       76652 : ldata_get_type(GEN ldata) { return mael3(ldata, 1, 1, 1); }
      39             : 
      40             : GEN
      41       74396 : ldata_get_an(GEN ldata) { return gel(ldata, 1); }
      42             : 
      43             : GEN
      44       59388 : ldata_get_dual(GEN ldata) { return gel(ldata, 2); }
      45             : 
      46             : long
      47        2813 : ldata_isreal(GEN ldata) { return isintzero(gel(ldata, 2)); }
      48             : 
      49             : GEN
      50      348495 : ldata_get_gammavec(GEN ldata) { return gel(ldata, 3); }
      51             : 
      52             : long
      53       26122 : ldata_get_degree(GEN ldata) { return lg(gel(ldata, 3))-1; }
      54             : 
      55             : GEN
      56      153825 : ldata_get_k(GEN ldata)
      57             : {
      58      153825 :   GEN w = gel(ldata,4);
      59      153825 :   if (typ(w) == t_VEC) w = gel(w,1);
      60      153825 :   return w;
      61             : }
      62             : 
      63             : /* a_n = O(n^{k1 + epsilon}) */
      64             : GEN
      65          98 : ldata_get_k1(GEN ldata)
      66             : {
      67          98 :   GEN w = gel(ldata,4);
      68          98 :   if (typ(w) == t_VEC) return gel(w,2);
      69             :   /* by default, assume that k1 = k-1 and even (k-1)/2 for entire functions */
      70          98 :   w = gaddgs(w,-1);
      71          98 :   return ldata_get_residue(ldata)? w: gmul2n(w, -1);
      72             : }
      73             : 
      74             : /* a_n = O(n^{k1 + epsilon}) */
      75             : static double
      76       91663 : ldata_get_k1_dbl(GEN ldata)
      77             : {
      78       91663 :   GEN w = gel(ldata,4);
      79             :   double k;
      80       91663 :   if (typ(w) == t_VEC) return gtodouble(gel(w,2));
      81             :   /* by default, assume that k1 = k-1 and even (k-1)/2 for entire functions */
      82       90228 :   k = gtodouble(w);
      83       90228 :   return ldata_get_residue(ldata)? k-1: (k-1)/2.;
      84             : }
      85             : 
      86             : GEN
      87      271868 : ldata_get_conductor(GEN ldata) { return gel(ldata, 5); }
      88             : 
      89             : GEN
      90      104930 : ldata_get_rootno(GEN ldata) { return gel(ldata, 6); }
      91             : 
      92             : GEN
      93      175164 : ldata_get_residue(GEN ldata) { return lg(ldata) == 7 ? NULL: gel(ldata, 7); }
      94             : 
      95             : long
      96      147266 : linit_get_type(GEN linit) { return mael(linit, 1, 1); }
      97             : 
      98             : GEN
      99      196969 : linit_get_ldata(GEN linit) { return gel(linit, 2); }
     100             : 
     101             : GEN
     102      249496 : linit_get_tech(GEN linit) { return gel(linit, 3); }
     103             : 
     104             : long
     105      299635 : is_linit(GEN data)
     106             : {
     107      184026 :   return lg(data) == 4 && typ(data) == t_VEC
     108      483661 :                        && typ(gel(data, 1)) == t_VECSMALL;
     109             : }
     110             : 
     111             : GEN
     112       31848 : lfun_get_step(GEN tech) { return gmael(tech, 2, 1);}
     113             : 
     114             : GEN
     115       31848 : lfun_get_pol(GEN tech) { return gmael(tech, 2, 2);}
     116             : 
     117             : GEN
     118        5151 : lfun_get_Residue(GEN tech) { return gmael(tech, 2, 3);}
     119             : 
     120             : GEN
     121       49928 : lfun_get_k2(GEN tech) { return gmael(tech, 3, 1);}
     122             : 
     123             : GEN
     124       19347 : lfun_get_w2(GEN tech) { return gmael(tech, 3, 2);}
     125             : 
     126             : GEN
     127       19347 : lfun_get_expot(GEN tech) { return gmael(tech, 3, 3);}
     128             : 
     129             : GEN
     130        9856 : lfun_get_factgammavec(GEN tech) { return gmael(tech, 3, 4); }
     131             : 
     132             : /* Handle complex Vga whose sum is real */
     133             : static GEN
     134      103486 : sumVga(GEN Vga) { return real_i(vecsum(Vga)); }
     135             : /* sum_i max (Im v[i],0) */
     136             : static double
     137       26980 : sumVgaimpos(GEN v)
     138             : {
     139       26980 :   double d = 0.;
     140       26980 :   long i, l = lg(v);
     141       74346 :   for (i = 1; i < l; i++)
     142             :   {
     143       47366 :     GEN c = imag_i(gel(v,i));
     144       47366 :     if (gsigne(c) > 0) d += gtodouble(c);
     145             :   }
     146       26980 :   return d;
     147             : }
     148             : 
     149             : static long
     150       42094 : vgaell(GEN Vga)
     151             : {
     152       42094 :   if (lg(Vga) == 3)
     153       30348 :   { GEN c = gsub(gel(Vga,1), gel(Vga,2)); return gequal1(c) || gequalm1(c); }
     154       11746 :   return 0;
     155             : }
     156             : int
     157       87181 : Vgaeasytheta(GEN Vga) { return lg(Vga)-1 == 1 || vgaell(Vga); }
     158             : /* return b(n) := a(n) * n^c, when Vgaeasytheta(Vga) is set */
     159             : static GEN
     160       18599 : antwist(GEN an, GEN Vga, long prec)
     161             : {
     162             :   long l, i;
     163       18599 :   GEN b, c = vecmin(Vga);
     164       18599 :   if (gequal0(c)) return an;
     165        4368 :   l = lg(an); b = cgetg(l, t_VEC);
     166        4368 :   if (gequal1(c))
     167             :   {
     168        3626 :     if (typ(an) == t_VECSMALL)
     169       17647 :       for (i = 1; i < l; i++) gel(b,i) = mulss(an[i], i);
     170             :     else
     171       41356 :       for (i = 1; i < l; i++) gel(b,i) = gmulgu(gel(an,i), i);
     172             :   }
     173             :   else
     174             :   {
     175         742 :     GEN v = vecpowug(l-1, c, prec);
     176         742 :     if (typ(an) == t_VECSMALL)
     177           0 :       for (i = 1; i < l; i++) gel(b,i) = gmulsg(an[i], gel(v,i));
     178             :     else
     179       33964 :       for (i = 1; i < l; i++) gel(b,i) = gmul(gel(an,i), gel(v,i));
     180             :   }
     181        4368 :   return b;
     182             : }
     183             : 
     184             : static GEN
     185        9618 : theta_dual(GEN theta, GEN bn)
     186             : {
     187        9618 :   if (typ(bn)==t_INT) return NULL;
     188             :   else
     189             :   {
     190          77 :     GEN thetad = shallowcopy(theta), ldata = linit_get_ldata(theta);
     191          77 :     GEN Vga = ldata_get_gammavec(ldata);
     192          77 :     GEN tech = shallowcopy(linit_get_tech(theta));
     193          77 :     GEN an = theta_get_an(tech);
     194          77 :     long prec = nbits2prec(theta_get_bitprec(tech));
     195          77 :     GEN vb = ldata_vecan(bn, lg(an)-1, prec);
     196          77 :     if (!theta_get_m(tech) && Vgaeasytheta(Vga)) vb = antwist(vb, Vga, prec);
     197          77 :     gel(tech,1) = vb;
     198          77 :     gel(thetad,3) = tech; return thetad;
     199             :   }
     200             : }
     201             : 
     202             : static GEN
     203       84370 : domain_get_dom(GEN domain)  { return gel(domain,1); }
     204             : static long
     205       25526 : domain_get_der(GEN domain)  { return mael2(domain, 2, 1); }
     206             : static long
     207       36824 : domain_get_bitprec(GEN domain)  { return mael2(domain, 2, 2); }
     208             : GEN
     209       84930 : lfun_get_domain(GEN tech) { return gel(tech,1); }
     210             : long
     211          91 : lfun_get_bitprec(GEN tech){ return domain_get_bitprec(lfun_get_domain(tech)); }
     212             : GEN
     213       59236 : lfun_get_dom(GEN tech) { return domain_get_dom(lfun_get_domain(tech)); }
     214             : 
     215             : GEN
     216        2575 : lfunprod_get_fact(GEN tech)  { return gel(tech, 2); }
     217             : 
     218             : GEN
     219       52454 : theta_get_an(GEN tdata)      { return gel(tdata, 1);}
     220             : GEN
     221        7889 : theta_get_K(GEN tdata)       { return gel(tdata, 2);}
     222             : GEN
     223        5991 : theta_get_R(GEN tdata)       { return gel(tdata, 3);}
     224             : long
     225       66012 : theta_get_bitprec(GEN tdata) { return itos(gel(tdata, 4));}
     226             : long
     227      102157 : theta_get_m(GEN tdata)       { return itos(gel(tdata, 5));}
     228             : GEN
     229       54015 : theta_get_tdom(GEN tdata)    { return gel(tdata, 6);}
     230             : GEN
     231       62884 : theta_get_isqrtN(GEN tdata)  { return gel(tdata, 7);}
     232             : 
     233             : /*******************************************************************/
     234             : /*  Helper functions related to Gamma products                     */
     235             : /*******************************************************************/
     236             : /* x != 0 */
     237             : static int
     238        6559 : serisscalar(GEN x)
     239             : {
     240             :   long i;
     241        6559 :   if (valser(x)) return 0;
     242        8218 :   for (i = lg(x)-1; i > 3; i--) if (!gequal0(gel(x,i))) return 0;
     243        6314 :   return 1;
     244             : }
     245             : 
     246             : /* return -itos(s) >= 0 if scalar s is (approximately) equal to a nonpositive
     247             :  * integer, and -1 otherwise */
     248             : static long
     249       20776 : isnegint(GEN s)
     250             : {
     251       20776 :   GEN r = ground(real_i(s));
     252       20776 :   if (signe(r) <= 0 && gequal(s, r)) return -itos(r);
     253       20650 :   return -1;
     254             : }
     255             : /* if s = a + O(x^n), a <= 0 integer, replace by a + b*x^n + O(x^(n+1)) */
     256             : static GEN
     257        6580 : serextendifnegint(GEN s, GEN b, long *ext)
     258             : {
     259        6580 :   if (!signe(s) || (serisscalar(s) && isnegint(gel(s,2)) >= 0))
     260             :   {
     261         112 :     long l = lg(s);
     262         112 :     GEN t = cgetg(l+1, t_SER);
     263         301 :     gel(t, l) = b; while (--l > 1) gel(t,l) = gel(s,l);
     264         112 :     if (gequal0(gel(t,2))) gel(t,2) = gen_0;
     265         112 :     t[1] = s[1]; s = normalizeser(t); *ext = 1;
     266             :   }
     267        6580 :   return s;
     268             : }
     269             : 
     270             : /* r/x + O(1), r != 0 */
     271             : static GEN
     272        5089 : serpole(GEN r)
     273             : {
     274        5089 :   GEN s = cgetg(3, t_SER);
     275        5089 :   s[1] = evalsigne(1)|evalvalser(-1)|evalvarn(0);
     276        5089 :   gel(s,2) = r; return s;
     277             : }
     278             : /* a0 +  a1 x + O(x^e), e >= 0 */
     279             : static GEN
     280        7511 : deg1ser_shallow(GEN a1, GEN a0, long v, long e)
     281        7511 : { return RgX_to_ser(deg1pol_shallow(a1, a0, v), e+2); }
     282             : 
     283             : /* pi^(-s/2) Gamma(s/2) */
     284             : static GEN
     285       10318 : gamma_R(GEN s, long *ext, long prec)
     286             : {
     287       10318 :   GEN s2 = gmul2n(s, -1);
     288             :   long ms;
     289             : 
     290       10318 :   if (typ(s) == t_SER)
     291        4816 :     s2 = serextendifnegint(s2, ghalf, ext);
     292        5502 :   else if ((ms = isnegint(s2)) >= 0)
     293             :   {
     294          35 :     GEN r = gmul(powPis(stoi(ms),prec), gdivsg(odd(ms)? -2: 2, mpfact(ms)));
     295          35 :     return serpole(r);
     296             :   }
     297       10283 :   return gdiv(ggamma(s2,prec), powPis(s2,prec));
     298             : }
     299             : /* gamma_R(s)gamma_R(s+1) = 2 (2pi)^(-s) Gamma(s) */
     300             : static GEN
     301       10724 : gamma_C(GEN s, long *ext, long prec)
     302             : {
     303             :   long ms;
     304       10724 :   if (typ(s) == t_SER)
     305        1764 :     s = serextendifnegint(s, gen_1, ext);
     306        8960 :   else if ((ms = isnegint(s)) >= 0)
     307             :   {
     308           0 :     GEN r = gmul(pow2Pis(stoi(ms),prec), gdivsg(odd(ms)? -2: 2, mpfact(ms)));
     309           0 :     return serpole(r);
     310             :   }
     311       10724 :   return gmul2n(gdiv(ggamma(s,prec), pow2Pis(s,prec)), 1);
     312             : }
     313             : 
     314             : static GEN
     315        1736 : gammafrac(GEN r, long d)
     316             : {
     317        1736 :   long i, l = labs(d) + 1, j = (d > 0)? 0: 2*d;
     318        1736 :   GEN T, v = cgetg(l, t_COL);
     319        6090 :   for (i = 1; i < l; i++, j += 2)
     320        4354 :     gel(v,i) = deg1pol_shallow(gen_1, gaddgs(r, j), 0);
     321        1736 :   T = RgV_prod(v); return d > 0? T: mkrfrac(gen_1, T);
     322             : }
     323             : 
     324             : /*
     325             : GR(s)=Pi^-(s/2)*gamma(s/2);
     326             : GC(s)=2*(2*Pi)^-s*gamma(s)
     327             : gdirect(F,s)=prod(i=1,#F,GR(s+F[i]))
     328             : gfact(F,s)=
     329             : { my([R,A,B]=gammafactor(F), [a,e]=A, [b,f]=B, p=poldegree(R));
     330             :   subst(R,x,s) * (2*Pi)^-p * prod(i=1,#a,GR(s+a[i])^e[i])
     331             :                            * prod(i=1,#b,GC(s+b[i])^f[i]); }
     332             : */
     333             : static GEN
     334       21791 : gammafactor(GEN Vga)
     335             : {
     336       21791 :   long i, r, c, l = lg(Vga);
     337       21791 :   GEN v, P, a, b, e, f, E, F = cgetg(l, t_VEC), R = gen_1;
     338       61278 :   for (i = 1; i < l; ++i)
     339             :   {
     340       39487 :     GEN a = gel(Vga,i), r = gmul2n(real_i(a), -1);
     341       39487 :     long q = itos(gfloor(r)); /* [Re a/2] */
     342       39487 :     r = gmul2n(gsubgs(r, q), 1);
     343       39487 :     gel(F,i) = gequal0(imag_i(a)) ? r : mkcomplex(r, imag_i(a)); /* 2{Re a/2} + I*(Im a) */
     344       39487 :     if (q) R = gmul(R, gammafrac(gel(F,i), q));
     345             :   }
     346       21791 :   F = vec_reduce(F, &E); l = lg(E);
     347       21791 :   v = cgetg(l, t_VEC);
     348       55314 :   for (i = 1; i < l; i++)
     349       33523 :       gel(v,i) = mkvec2(gsub(gel(F,i),gfloor(real_i(gel(F,i)))), stoi(E[i]));
     350       21791 :   gen_sort_inplace(v, (void*)cmp_universal, cmp_nodata, &P);
     351       21791 :   a = cgetg(l, t_VEC); e = cgetg(l, t_VECSMALL);
     352       21791 :   b = cgetg(l, t_VEC); f = cgetg(l, t_VECSMALL);
     353       44548 :   for (i = r = c = 1; i < l;)
     354       22757 :     if (i==l-1 || cmp_universal(gel(v,i), gel(v,i+1)))
     355       11991 :     { gel(a, r) = gel(F, P[i]); e[r++] = E[P[i]]; i++; }
     356             :     else
     357       10766 :     { gel(b, c) = gel(F, P[i]); f[c++] = E[P[i]]; i+=2; }
     358       21791 :   setlg(a, r); setlg(e, r);
     359       21791 :   setlg(b, c); setlg(f, c); return mkvec3(R, mkvec2(a,e), mkvec2(b,f));
     360             : }
     361             : 
     362             : static GEN
     363        3654 : polgammaeval(GEN F, GEN s)
     364             : {
     365        3654 :   GEN r = poleval(F, s);
     366        3654 :   if (typ(s) != t_SER && gequal0(r))
     367             :   { /* here typ(F) = t_POL */
     368             :     long e;
     369           7 :     for (e = 1;; e++)
     370             :     {
     371           7 :       F = RgX_deriv(F); r = poleval(F,s);
     372           7 :       if (!gequal0(r)) break;
     373             :     }
     374           7 :     if (e > 1) r = gdiv(r, mpfact(e));
     375           7 :     r = serpole(r); setvalser(r, e);
     376             :   }
     377        3654 :   return r;
     378             : }
     379             : static long
     380        1799 : rfrac_degree(GEN R)
     381             : {
     382        1799 :   GEN a = gel(R,1), b = gel(R,2);
     383        1799 :   return ((typ(a) == t_POL)? degpol(a): 0) - degpol(b);
     384             : }
     385             : static GEN
     386       20048 : fracgammaeval(GEN F, GEN s, long prec)
     387             : {
     388       20048 :   GEN R = gel(F,1);
     389             :   long d;
     390       20048 :   switch(typ(R))
     391             :   {
     392          56 :     case t_POL:
     393          56 :       d = degpol(R);
     394          56 :       R = polgammaeval(R, s); break;
     395        1799 :     case t_RFRAC:
     396        1799 :       d = rfrac_degree(R);
     397        1799 :       R = gdiv(polgammaeval(gel(R,1), s), polgammaeval(gel(R,2), s)); break;
     398       18193 :     default: return R;
     399             :   }
     400        1855 :   return gmul(R, powrs(Pi2n(1,prec), -d));
     401             : }
     402             : 
     403             : static GEN
     404       20048 : gammafactproduct(GEN F, GEN s, long *ext, long prec)
     405             : {
     406       20048 :   pari_sp av = avma;
     407       20048 :   GEN R = gel(F,2), Rw = gel(R,1), Re = gel(R,2);
     408       20048 :   GEN C = gel(F,3), Cw = gel(C,1), Ce = gel(C,2), z = fracgammaeval(F,s,prec);
     409       20048 :   long i, lR = lg(Rw), lC = lg(Cw);
     410       20048 :   *ext = 0;
     411       30366 :   for (i = 1; i < lR; i++)
     412       10318 :     z = gmul(z, gpowgs(gamma_R(gadd(s,gel(Rw, i)), ext, prec), Re[i]));
     413       30772 :   for (i = 1; i < lC; i++)
     414       10724 :     z = gmul(z, gpowgs(gamma_C(gadd(s,gel(Cw, i)), ext, prec), Ce[i]));
     415       20048 :   return gc_upto(av, z);
     416             : }
     417             : 
     418             : static int
     419        5166 : gammaordinary(GEN Vga, GEN s)
     420             : {
     421        5166 :   long i, d = lg(Vga)-1;
     422       13839 :   for (i = 1; i <= d; i++)
     423             :   {
     424        8764 :     GEN z = gadd(s, gel(Vga,i));
     425             :     long e;
     426        8764 :     if (gexpo(imag_i(z)) < -10)
     427             :     {
     428        8764 :       z = real_i(z);
     429        8764 :       if (gsigne(z) <= 0) { (void)grndtoi(z, &e); if (e < -10) return 0; }
     430             :     }
     431             :   }
     432        5075 :   return 1;
     433             : }
     434             : 
     435             : /* Exponent A of t in asymptotic expansion; K(t) ~ C t^A exp(-pi d t^(2/d)).
     436             :  * suma = vecsum(Vga)*/
     437             : static double
     438       91656 : gammavec_expo(long d, double suma) { return (1 - d + suma) / d; }
     439             : 
     440             : /*******************************************************************/
     441             : /*       First part: computations only involving Theta(t)          */
     442             : /*******************************************************************/
     443             : 
     444             : static void
     445      139761 : get_cone(GEN t, double *r, double *a)
     446             : {
     447      139761 :   const long prec = LOWDEFAULTPREC;
     448      139761 :   if (typ(t) == t_COMPLEX)
     449             :   {
     450       21532 :     t  = gprec_w(t, prec);
     451       21532 :     *r = gtodouble(gabs(t, prec));
     452       21532 :     *a = fabs(gtodouble(garg(t, prec)));
     453             :   }
     454             :   else
     455             :   {
     456      118229 :     *r = fabs(gtodouble(t));
     457      118229 :     *a = 0.;
     458             :   }
     459      139761 :   if (!*r && !*a) pari_err_DOMAIN("lfunthetainit","t","=",gen_0,t);
     460      139754 : }
     461             : /* slightly larger cone than necessary, to avoid round-off problems */
     462             : static void
     463       85746 : get_cone_fuzz(GEN t, double *r, double *a)
     464       85746 : { get_cone(t, r, a); *r -= 1e-10; if (*a) *a += 1e-10; }
     465             : 
     466             : /* Initialization m-th Theta derivative. tdom is either
     467             :  * - [rho,alpha]: assume |t| >= rho and |arg(t)| <= alpha
     468             :  * - a positive real scalar: assume t real, t >= tdom;
     469             :  * - a complex number t: compute at t;
     470             :  * N is the conductor (either the true one from ldata or a guess from
     471             :  * lfunconductor) */
     472             : long
     473       64683 : lfunthetacost(GEN ldata, GEN tdom, long m, long bit, long *extrabit)
     474             : {
     475       64683 :   pari_sp av = avma;
     476       64683 :   GEN Vga = ldata_get_gammavec(ldata);
     477       64683 :   long d = lg(Vga)-1;
     478       64683 :   double k1 = maxdd(ldata_get_k1_dbl(ldata), 0.);
     479       64683 :   double c = d/2., a, A, B, logC, al, rho, T;
     480       64683 :   double N = gtodouble(ldata_get_conductor(ldata));
     481             : 
     482       64683 :   if (extrabit) *extrabit = 0;
     483       64683 :   if (!N) pari_err_TYPE("lfunthetaneed [missing conductor]", ldata);
     484       64683 :   if (typ(tdom) == t_VEC && lg(tdom) == 3)
     485             :   {
     486           7 :     rho= gtodouble(gel(tdom,1));
     487           7 :     al = gtodouble(gel(tdom,2));
     488             :   }
     489             :   else
     490       64676 :     get_cone_fuzz(tdom, &rho, &al);
     491       64676 :   A = gammavec_expo(d, gtodouble(sumVga(Vga))); set_avma(av);
     492       64676 :   a = (A+k1+1) + (m-1)/c;
     493       64676 :   if (fabs(a) < 1e-10) a = 0.;
     494       64676 :   logC = c*M_LN2 - log(c)/2;
     495             :   /* +1: fudge factor */
     496       64676 :   B = M_LN2*bit+logC+m*log(2*M_PI) + 1 + (k1+1)*log(N)/2 - (k1+m+1)*log(rho);
     497       64676 :   if (al)
     498             :   { /* t = rho e^(i*al), T^(1/c) = Re(t^(1/c)) > 0, T = rho cos^c(al/c) */
     499       10773 :     double z = cos(al/c);
     500       10773 :     if (z <= 0)
     501           7 :       pari_err_DOMAIN("lfunthetaneed", "arg t", ">", dbltor(c*M_PI/2), tdom);
     502       10766 :     T = (d == 2 && typ(tdom) != t_VEC)? gtodouble(real_i(tdom)): rho*pow(z,c);
     503       10766 :     B -= log(z) * (c * (k1+A+1) + m);
     504             :   }
     505             :   else
     506       53903 :     T = rho;
     507       64669 :   if (B <= 0) return 0;
     508       64669 :   A = floor(0.9 + dblcoro526(a,c,B) / T * sqrt(N));
     509       64669 :   if (dblexpo(A) >= BITS_IN_LONG-1) pari_err_OVERFLOW("lfunthetacost");
     510       64662 :   if (extrabit && a * log2(A) > 16) *extrabit = a * log2(A);
     511       64662 :   return (long)A;
     512             : }
     513             : long
     514          21 : lfunthetacost0(GEN L, GEN tdom, long m, long bitprec)
     515             : {
     516             :   long n;
     517          21 :   if (is_linit(L) && linit_get_type(L)==t_LDESC_THETA)
     518           7 :   {
     519           7 :     GEN tech = linit_get_tech(L);
     520           7 :     n = lg(theta_get_an(tech))-1;
     521             :   }
     522             :   else
     523             :   {
     524          14 :     pari_sp av = avma;
     525          14 :     GEN ldata = lfunmisc_to_ldata_shallow(L);
     526          14 :     n = lfunthetacost(ldata, tdom? tdom: gen_1, m, bitprec, NULL);
     527           7 :     set_avma(av);
     528             :   }
     529          14 :   return n;
     530             : }
     531             : 
     532             : static long
     533        6349 : fracgammadegree(GEN FVga)
     534        6349 : { GEN F = gel(FVga,1); return (typ(F)==t_RFRAC)? degpol(gel(F,2)): 0; }
     535             : 
     536             : /* Poles of a L-function can be represented in the following ways:
     537             :  * 1) Nothing (ldata has only 6 components, ldata_get_residue = NULL).
     538             :  * 2) a complex number (single pole at s = k with given residue, unknown if 0).
     539             :  * 3) A vector (possibly empty) of 2-component vectors [a, ra], where a is the
     540             :  * pole, ra a t_SER: its Taylor expansion at a. A t_VEC encodes the polar
     541             :  * part of L, a t_COL, the polar part of Lambda */
     542             : 
     543             : /* 'a' a complex number (pole), 'r' the polar part of L at 'a';
     544             :  * return 'R' the polar part of Lambda at 'a' */
     545             : static GEN
     546        4704 : rtoR(GEN a, GEN r, GEN FVga, GEN N, long prec)
     547             : {
     548        4704 :   long v = lg(r)-2, d = fracgammadegree(FVga), ext;
     549        4704 :   GEN Na, as = deg1ser_shallow(gen_1, a, varn(r), v);
     550        4704 :   Na = gpow(N, gdivgu(as, 2), prec);
     551             :   /* make up for a possible loss of accuracy */
     552        4704 :   if (d) as = deg1ser_shallow(gen_1, a, varn(r), v + d);
     553        4704 :   return gmul(gmul(r, Na), gammafactproduct(FVga, as, &ext, prec));
     554             : }
     555             : 
     556             : /* assume r in normalized form: t_VEC of pairs [be,re] */
     557             : GEN
     558        4473 : lfunrtopoles(GEN r)
     559             : {
     560        4473 :   long j, l = lg(r);
     561        4473 :   GEN v = cgetg(l, t_VEC);
     562        9177 :   for (j = 1; j < l; j++)
     563             :   {
     564        4704 :     GEN rj = gel(r,j), a = gel(rj,1);
     565        4704 :     gel(v,j) = a;
     566             :   }
     567        4473 :   gen_sort_inplace(v, (void*)&cmp_universal, cmp_nodata, NULL);
     568        4473 :   return v;
     569             : }
     570             : 
     571             : /* r / x + O(1) */
     572             : static GEN
     573        5194 : simple_pole(GEN r)
     574        5194 : { return isintzero(r)? gen_0: serpole(r); }
     575             : static GEN
     576        6048 : normalize_simple_pole(GEN r, GEN k)
     577             : {
     578        6048 :   long tx = typ(r);
     579        6048 :   if (is_vec_t(tx)) return r;
     580        5194 :   if (!is_scalar_t(tx)) pari_err_TYPE("lfunrootres [poles]", r);
     581        5194 :   return mkvec(mkvec2(k, simple_pole(r)));
     582             : }
     583             : /* normalize the description of a polar part */
     584             : static GEN
     585        5348 : normalizepoles(GEN r, GEN k)
     586             : {
     587             :   long iv, j, l;
     588             :   GEN v;
     589        5348 :   if (!is_vec_t(typ(r))) return normalize_simple_pole(r, k);
     590        2233 :   v = cgetg_copy(r, &l);
     591        5593 :   for (j = iv = 1; j < l; j++)
     592             :   {
     593        3360 :     GEN rj = gel(r,j), a = gel(rj,1), ra = gel(rj,2);
     594        3360 :     if (!is_scalar_t(typ(a)) || typ(ra) != t_SER)
     595           0 :       pari_err_TYPE("lfunrootres [poles]",r);
     596        3360 :     if (valser(ra) >= 0) continue;
     597        3346 :     gel(v,iv++) = rj;
     598             :   }
     599        2233 :   setlg(v, iv); return v;
     600             : }
     601             : static int
     602        8694 : residues_known(GEN r)
     603             : {
     604        8694 :   long i, l = lg(r);
     605        8694 :   if (isintzero(r)) return 0;
     606        8365 :   if (!is_vec_t(typ(r))) return 1;
     607       11088 :   for (i = 1; i < l; i++)
     608             :   {
     609        6748 :     GEN ri = gel(r,i);
     610        6748 :     if (!is_vec_t(typ(ri)) || lg(ri)!=3)
     611           0 :       pari_err_TYPE("lfunrootres [poles]",r);
     612        6748 :     if (isintzero(gel(ri, 2))) return 0;
     613             :   }
     614        4340 :   return 1;
     615             : }
     616             : 
     617             : /* Compute R's from r's (r = Taylor devts of L(s), R of Lambda(s)).
     618             :  * 'r/eno' passed to override the one from ldata  */
     619             : static GEN
     620       23674 : lfunrtoR_i(GEN ldata, GEN r, GEN eno, long prec)
     621             : {
     622       23674 :   GEN Vga = ldata_get_gammavec(ldata), N = ldata_get_conductor(ldata);
     623             :   GEN R, vr, FVga;
     624       23674 :   pari_sp av = avma;
     625             :   long lr, j, jR;
     626       23674 :   GEN k = ldata_get_k(ldata);
     627             : 
     628       23674 :   if (!r || isintzero(eno) || !residues_known(r))
     629       18326 :     return gen_0;
     630        5348 :   r = normalizepoles(r, k);
     631        5348 :   if (typ(r) == t_COL) return gc_GEN(av, r);
     632        4473 :   if (typ(ldata_get_dual(ldata)) != t_INT)
     633           0 :     pari_err(e_MISC,"please give the Taylor expansion of Lambda");
     634        4473 :   vr = lfunrtopoles(r); lr = lg(vr);
     635        4473 :   FVga = gammafactor(Vga);
     636        4473 :   R = cgetg(2*lr, t_COL);
     637        9177 :   for (j = jR = 1; j < lr; j++)
     638             :   {
     639        4704 :     GEN rj = gel(r,j), a = gel(rj,1), ra = gel(rj,2);
     640        4704 :     GEN Ra = rtoR(a, ra, FVga, N, prec);
     641        4704 :     GEN b = gsub(k, conj_i(a));
     642        4704 :     if (lg(Ra)-2 < -valser(Ra))
     643           0 :       pari_err(e_MISC,
     644             :         "please give more terms in L function's Taylor expansion at %Ps", a);
     645        4704 :     gel(R,jR++) = mkvec2(a, Ra);
     646        4704 :     if (!tablesearch(vr, b, (int (*)(GEN,GEN))&cmp_universal))
     647             :     {
     648        4515 :       GEN mX = gneg(pol_x(varn(Ra)));
     649        4515 :       GEN Rb = gmul(eno, gsubst(conj_i(Ra), varn(Ra), mX));
     650        4515 :       gel(R,jR++) = mkvec2(b, Rb);
     651             :     }
     652             :   }
     653        4473 :   setlg(R, jR); return gc_GEN(av, R);
     654             : }
     655             : static GEN
     656       23198 : lfunrtoR_eno(GEN ldata, GEN eno, long prec)
     657       23198 : { return lfunrtoR_i(ldata, ldata_get_residue(ldata), eno, prec); }
     658             : static GEN
     659       21077 : lfunrtoR(GEN ldata, long prec)
     660       21077 : { return lfunrtoR_eno(ldata, ldata_get_rootno(ldata), prec); }
     661             : 
     662             : static long
     663       21077 : prec_fix(long prec)
     664             : {
     665             : #ifndef LONG_IS_64BIT
     666             :   /* make sure that default accuracy is the same on 32/64bit */
     667        3011 :   if (odd(prec)) prec += EXTRAPREC64;
     668             : #endif
     669       21077 :   return prec;
     670             : }
     671             : 
     672             : /* thetainit using {an: n <= L}; if (m = 0 && easytheta), an2 is an * n^al */
     673             : static GEN
     674       21077 : lfunthetainit0(GEN ldata, GEN tdom, GEN an2, long m,
     675             :     long bitprec, long extrabit)
     676             : {
     677       21077 :   long prec = nbits2prec(bitprec);
     678       21077 :   GEN tech, N = ldata_get_conductor(ldata);
     679       21077 :   GEN K = gammamellininvinit(ldata, m, bitprec + extrabit);
     680       21077 :   GEN R = lfunrtoR(ldata, prec);
     681       21077 :   if (!tdom) tdom = gen_1;
     682       21077 :   if (typ(tdom) != t_VEC)
     683             :   {
     684             :     double r, a;
     685       21070 :     get_cone_fuzz(tdom, &r, &a);
     686       21070 :     tdom = mkvec2(dbltor(r), a? dbltor(a): gen_0);
     687             :   }
     688       21077 :   prec += maxss(EXTRAPREC64, nbits2extraprec(extrabit));
     689       21077 :   tech = mkvecn(7, an2,K,R, stoi(bitprec), stoi(m), tdom,
     690             :                    gsqrt(ginv(N), prec_fix(prec)));
     691       21077 :   return mkvec3(mkvecsmall(t_LDESC_THETA), ldata, tech);
     692             : }
     693             : 
     694             : /* tdom: 1) positive real number r, t real, t >= r; or
     695             :  *       2) [r,a], describing the cone |t| >= r, |arg(t)| <= a */
     696             : static GEN
     697       10367 : lfunthetainit_i(GEN data, GEN tdom, long m, long bit)
     698             : {
     699       10367 :   GEN ldata = lfunmisc_to_ldata_shallow(data);
     700       10367 :   long extrabit, b = 32, L = lfunthetacost(ldata, tdom, m, bit, &extrabit);
     701       10353 :   long prec = nbits2prec(bit + extrabit);
     702       10353 :   GEN ldatan = ldata_newprec(ldata, prec);
     703       10353 :   GEN an = ldata_get_an(ldatan), Vga = ldata_get_gammavec(ldatan);
     704       10353 :   an = ldata_vecan(an, L, prec);
     705       10353 :   if (m == 0 && Vgaeasytheta(Vga)) an = antwist(an, Vga, prec);
     706       10353 :   if (typ(an) != t_VECSMALL) b = maxss(b, gexpo(an));
     707       10353 :   return lfunthetainit0(ldatan, tdom, an, m, bit, b);
     708             : }
     709             : 
     710             : GEN
     711         336 : lfunthetainit(GEN ldata, GEN tdom, long m, long bitprec)
     712             : {
     713         336 :   pari_sp av = avma;
     714         336 :   GEN S = lfunthetainit_i(ldata, tdom? tdom: gen_1, m, bitprec);
     715         336 :   return gc_GEN(av, S);
     716             : }
     717             : 
     718             : GEN
     719        2450 : lfunan(GEN ldata, long L, long prec)
     720             : {
     721        2450 :   pari_sp av = avma;
     722             :   GEN an ;
     723        2450 :   ldata = ldata_newprec(lfunmisc_to_ldata_shallow(ldata), prec);
     724        2450 :   an = gc_GEN(av, ldata_vecan(ldata_get_an(ldata), L, prec));
     725        2394 :   if (typ(an) != t_VEC) an = vecsmall_to_vec_inplace(an);
     726        2394 :   return an;
     727             : }
     728             : 
     729             : static GEN
     730       15897 : mulrealvec(GEN x, GEN y)
     731             : {
     732       15897 :   if (is_vec_t(typ(x)) && is_vec_t(typ(y)))
     733          84 :     pari_APPLY_same(mulreal(gel(x,i),gel(y,i)))
     734             :   else
     735       15869 :     return mulreal(x,y);
     736             : }
     737             : static GEN
     738       31449 : gmulvec(GEN x, GEN y)
     739             : {
     740       31449 :   if (is_vec_t(typ(x)) && is_vec_t(typ(y)))
     741        2702 :     pari_APPLY_same(gmul(gel(x,i),gel(y,i)))
     742             :   else
     743       30784 :     return gmul(x,y);
     744             : }
     745             : static GEN
     746        9597 : gdivvec(GEN x, GEN y)
     747             : {
     748        9597 :   if (is_vec_t(typ(x)) && is_vec_t(typ(y)))
     749        2247 :     pari_APPLY_same(gdiv(gel(x,i),gel(y,i)))
     750             :   else
     751        9023 :     return gdiv(x,y);
     752             : }
     753             : 
     754             : static GEN
     755        3556 : gsubvec(GEN x, GEN y)
     756             : {
     757        3556 :   if (is_vec_t(typ(x)) && !is_vec_t(typ(y)))
     758           0 :     pari_APPLY_same(gsub(gel(x,i),y))
     759             :   else
     760        3556 :     return gsub(x,y);
     761             : }
     762             : 
     763             : /* return [1^(2/d), 2^(2/d),...,lim^(2/d)] */
     764             : static GEN
     765        7889 : mkvroots(long d, long lim, long prec)
     766             : {
     767        7889 :   if (d <= 4)
     768             :   {
     769        7539 :     GEN v = cgetg(lim+1,t_VEC);
     770             :     long n;
     771        7539 :     switch(d)
     772             :     {
     773        2800 :       case 1:
     774       52936 :         for (n=1; n <= lim; n++) gel(v,n) = sqru(n);
     775        2800 :         return v;
     776        1281 :       case 2:
     777      228529 :         for (n=1; n <= lim; n++) gel(v,n) = utoipos(n);
     778        1281 :         return v;
     779        1990 :       case 4:
     780     6117579 :         for (n=1; n <= lim; n++) gel(v,n) = sqrtr(utor(n, prec));
     781        1990 :         return v;
     782             :     }
     783             :   }
     784        1818 :   return vecpowug(lim, gdivgu(gen_2,d), prec);
     785             : }
     786             : 
     787             : GEN
     788       62114 : lfunthetacheckinit(GEN data, GEN t, long m, long bitprec)
     789             : {
     790       62114 :   if (is_linit(data) && linit_get_type(data)==t_LDESC_THETA)
     791             :   {
     792       54015 :     GEN tdom, thetainit = linit_get_tech(data);
     793       54015 :     long bitprecnew = theta_get_bitprec(thetainit);
     794       54015 :     long m0 = theta_get_m(thetainit);
     795             :     double r, al, rt, alt;
     796       54015 :     if (m0 != m)
     797           0 :       pari_err_DOMAIN("lfuntheta","derivative order","!=", stoi(m),stoi(m0));
     798       54015 :     if (bitprec > bitprecnew) goto INIT;
     799       54015 :     get_cone(t, &rt, &alt);
     800       54015 :     tdom = theta_get_tdom(thetainit);
     801       54015 :     r = gtodouble(gel(tdom,1));
     802       54015 :     al= gtodouble(gel(tdom,2)); if (rt >= r && alt <= al) return data;
     803             :   }
     804        8099 : INIT:
     805        9940 :   return lfunthetainit_i(data, t, m, bitprec);
     806             : }
     807             : 
     808             : static GEN
     809    14712655 : get_an(GEN an, long n)
     810             : {
     811    14712655 :   if (typ(an) == t_VECSMALL) { long a = an[n]; if (a) return stoi(a); }
     812    14712655 :   else { GEN a = gel(an,n); if (a && !gequal0(a)) return a; }
     813    12685598 :   return NULL;
     814             : }
     815             : /* x * an[n] */
     816             : static GEN
     817    12877655 : mul_an(GEN an, long n, GEN x)
     818             : {
     819    12877655 :   if (typ(an) == t_VECSMALL) { long a = an[n]; if (a) return gmulsg(a,x); }
     820     7455528 :   else { GEN a = gel(an,n); if (a && !gequal0(a)) return gmul(a,x); }
     821     2440507 :   return NULL;
     822             : }
     823             : /* 2*t^a * x **/
     824             : static GEN
     825      334909 : mulT(GEN t, GEN a, GEN x, long prec)
     826             : {
     827      334909 :   if (gequal0(a)) return gmul2n(x,1);
     828       32647 :   return gmul(x, gmul2n(gequal1(a)? t: gpow(t,a,prec), 1));
     829             : }
     830             : 
     831             : static GEN
     832    35381321 : vecan_cmul(void *E, GEN P, long a, GEN x)
     833             : {
     834             :   (void)E;
     835    35381321 :   if (typ(P) == t_VECSMALL)
     836    24516803 :     return (a==0 || !P[a])? NULL: gmulsg(P[a], x);
     837             :   else
     838    10864518 :     return (a==0 || !gel(P,a))? NULL: gmul(gel(P,a), x);
     839             : }
     840             : /* d=2, 2 sum_{n <= N} a(n) (n t)^al q^n, q = exp(-2pi t),
     841             :  * an2[n] = a(n) * n^al */
     842             : static GEN
     843      294984 : theta2_i(GEN an2, long N, GEN t, GEN al, long prec)
     844             : {
     845      294984 :   GEN S, q, pi2 = Pi2n(1,prec);
     846      294978 :   const struct bb_algebra *alg = get_Rg_algebra();
     847      294978 :   setsigne(pi2,-1); q = gexp(gmul(pi2, t), prec);
     848             :   /* Brent-Kung in case the a_n are small integers */
     849      294979 :   S = gen_bkeval(an2, N, q, 1, NULL, alg, vecan_cmul);
     850      294981 :   return mulT(t, al, S, prec);
     851             : }
     852             : static GEN
     853      286847 : theta2(GEN an2, long N, GEN t, GEN al, long prec)
     854             : {
     855      286847 :   pari_sp av = avma;
     856      286847 :   return gc_upto(av, theta2_i(an2, N, t, al, prec));
     857             : }
     858             : 
     859             : /* d=1, 2 sum_{n <= N} a_n (n t)^al q^(n^2), q = exp(-pi t^2),
     860             :  * an2[n] is a_n n^al */
     861             : static GEN
     862       39928 : theta1(GEN an2, long N, GEN t, GEN al, long prec)
     863             : {
     864       39928 :   GEN q = gexp(gmul(negr(mppi(prec)), gsqr(t)), prec);
     865       39928 :   GEN vexp = gsqrpowers(q, N), S = gen_0;
     866       39928 :   pari_sp av = avma;
     867             :   long n;
     868     6783691 :   for (n = 1; n <= N; n++)
     869             :   {
     870     6743763 :     GEN c = mul_an(an2, n, gel(vexp,n));
     871     6743763 :     if (c)
     872             :     {
     873     5667353 :       S = gadd(S, c);
     874     5667353 :       if (gc_needed(av, 3)) S = gc_upto(av, S);
     875             :     }
     876             :   }
     877       39928 :   return mulT(t, al, S, prec);
     878             : }
     879             : 
     880             : /* If m > 0, compute m-th derivative of theta(t) = theta0(t/sqrt(N))
     881             :  * with absolute error 2^-bitprec; theta(t)=\sum_{n\ge1}a(n)K(nt/N^(1/2)) */
     882             : GEN
     883       52167 : lfuntheta(GEN data, GEN t, long m, long bitprec)
     884             : {
     885       52167 :   pari_sp ltop = avma;
     886             :   long limt, d;
     887             :   GEN isqN, vecan, Vga, ldata, theta, thetainit, S;
     888             :   long n, prec;
     889             : 
     890       52167 :   theta = lfunthetacheckinit(data, t, m, bitprec);
     891       52160 :   ldata = linit_get_ldata(theta);
     892       52160 :   thetainit = linit_get_tech(theta);
     893       52160 :   vecan = theta_get_an(thetainit);
     894       52160 :   isqN = theta_get_isqrtN(thetainit);
     895       52160 :   prec = maxss(realprec(isqN), nbits2prec(bitprec));
     896       52160 :   t = gprec_w(t, prec);
     897       52160 :   limt = lg(vecan)-1;
     898       52160 :   if (theta == data)
     899       47967 :     limt = minss(limt, lfunthetacost(ldata, t, m, bitprec, NULL));
     900       52160 :   if (!limt)
     901             :   {
     902          14 :     set_avma(ltop); S = real_0_bit(-bitprec);
     903          14 :     if (!is_real_t(typ(t)) || !ldata_isreal(ldata))
     904           7 :       S = gc_GEN(ltop, mkcomplex(S,S));
     905          14 :     return S;
     906             :   }
     907       52146 :   t = gmul(t, isqN);
     908       52146 :   Vga = ldata_get_gammavec(ldata);
     909       52146 :   d = lg(Vga)-1;
     910       52146 :   if (m == 0 && Vgaeasytheta(Vga))
     911             :   {
     912       48065 :     if (theta_get_m(thetainit) > 0) vecan = antwist(vecan, Vga, prec);
     913       48065 :     if (d == 1) S = theta1(vecan, limt, t, gel(Vga,1), prec);
     914        8137 :     else        S = theta2_i(vecan, limt, t, vecmin(Vga), prec);
     915             :   }
     916             :   else
     917             :   {
     918        4081 :     GEN K = theta_get_K(thetainit);
     919        4081 :     GEN vroots = mkvroots(d, limt, prec);
     920             :     pari_sp av;
     921        4081 :     t = gpow(t, gdivgu(gen_2,d), prec);
     922        4081 :     S = gen_0; av = avma;
     923    14716736 :     for (n = 1; n <= limt; ++n)
     924             :     {
     925    14712655 :       GEN nt, an = get_an(vecan, n);
     926    14712655 :       if (!an) continue;
     927     2027057 :       nt = gmul(gel(vroots,n), t);
     928     2027057 :       if (m) an = gmul(an, powuu(n, m));
     929     2027057 :       S = gadd(S, gmul(an, gammamellininvrt(K, nt, bitprec)));
     930     2027057 :       if ((n & 0x1ff) == 0) S = gc_upto(av, S);
     931             :     }
     932        4081 :     if (m) S = gmul(S, gpowgs(isqN, m));
     933             :   }
     934       52146 :   return gc_upto(ltop, S);
     935             : }
     936             : 
     937             : /*******************************************************************/
     938             : /* Second part: Computation of L-Functions.                        */
     939             : /*******************************************************************/
     940             : 
     941             : struct lfunp {
     942             :   long precmax, Dmax, D, M, m0, nmax, d, vgaell;
     943             :   double k1, dc, dw, dh, MAXs, sub;
     944             :   GEN L, an, bn;
     945             : };
     946             : 
     947             : static void
     948       26980 : lfunp_set(GEN ldata, long der, long bitprec, struct lfunp *S)
     949             : {
     950       26980 :   const long derprec = (der > 1)? dbllog2(mpfact(der)): 0; /* log2(der!) */
     951             :   GEN Vga, N, L, k;
     952             :   long k1, d, m, M, flag, nmax;
     953             :   double a, A, E, hd, Ep, d2, suma, maxs, mins, sub, B0,B1;
     954             :   double logN2, logC, Lestimate, Mestimate;
     955             : 
     956       26980 :   Vga = ldata_get_gammavec(ldata);
     957       26980 :   S->d = d = lg(Vga)-1; d2 = d/2.;
     958             : 
     959       26980 :   suma = gtodouble(sumVga(Vga));
     960       26980 :   k = ldata_get_k(ldata);
     961       26980 :   N = ldata_get_conductor(ldata);
     962       26980 :   logN2 = log(gtodouble(N)) / 2;
     963       26980 :   maxs = S->dc + S->dw;
     964       26980 :   mins = S->dc - S->dw;
     965       26980 :   S->MAXs = maxdd(maxs, gtodouble(k)-mins);
     966             : 
     967             :   /* we compute Lambda^(der)(s) / der!; need to compensate for L^(der)(s)
     968             :    * ln |gamma(s)| ~ -(pi/4) \sum_i |Im(s + a_i)|; max with 1: fudge factor */
     969       26980 :   a = (M_PI/(4*M_LN2))*(d*S->dh + sumVgaimpos(Vga));
     970       26980 :   S->D = (long)ceil(bitprec + derprec + maxdd(a, 1));
     971       26980 :   E = M_LN2*S->D; /* D:= required absolute bitprec */
     972             : 
     973       26980 :   Ep = E + maxdd(M_PI * S->dh * d2, (d*S->MAXs + suma - 1) * log(E));
     974       26980 :   hd = d2*M_PI*M_PI / Ep;
     975       26980 :   S->m0 = (long)ceil(M_LN2/hd);
     976       26980 :   hd = M_LN2/S->m0;
     977             : 
     978       26980 :   logC = d2*M_LN2 - log(d2)/2;
     979       26980 :   k1 = maxdd(ldata_get_k1_dbl(ldata), 0.);
     980       26980 :   S->k1 = k1; /* assume |a_n| << n^k1 with small implied constant */
     981       26980 :   A = gammavec_expo(d, suma);
     982             : 
     983       26980 :   sub = 0.;
     984       26980 :   if (mins > 1)
     985             :   {
     986        5166 :     GEN sig = dbltor(mins);
     987        5166 :     sub += logN2*mins;
     988        5166 :     if (gammaordinary(Vga, sig))
     989             :     {
     990             :       long ext;
     991        5075 :       GEN gas = gammafactproduct(gammafactor(Vga), sig, &ext, LOWDEFAULTPREC);
     992        5075 :       if (typ(gas) != t_SER)
     993             :       {
     994        5075 :         double dg = dbllog2(gas);
     995        5075 :         if (dg > 0) sub += dg * M_LN2;
     996             :       }
     997             :     }
     998             :   }
     999       26980 :   S->sub = sub;
    1000       26980 :   M = 1000;
    1001       26980 :   L = cgetg(M+2, t_VECSMALL);
    1002       26980 :   a = S->k1 + A;
    1003             : 
    1004       26980 :   B0 = 5 + E - S->sub + logC + S->k1*logN2; /* 5 extra bits */
    1005       26980 :   B1 = hd * (S->MAXs - S->k1);
    1006       26980 :   Lestimate = dblcoro526(a + S->MAXs - 2./d, d/2.,
    1007       26980 :     E - S->sub + logC - log(2*M_PI*hd) + S->MAXs*logN2);
    1008       26980 :   Mestimate = ((Lestimate > 0? log(Lestimate): 0) + logN2) / hd;
    1009       26980 :   nmax = 0;
    1010       26980 :   flag = 0;
    1011       26980 :   for (m = 0;; m++)
    1012     2414202 :   {
    1013     2441182 :     double x, H = logN2 - m*hd, B = B0 + m*B1;
    1014             :     long n;
    1015     2441182 :     x = dblcoro526(a, d/2., B);
    1016     2441182 :     n = floor(x*exp(H));
    1017     2441182 :     if (n > nmax) nmax = n;
    1018     2441182 :     if (m > M) { M *= 2; L = vecsmall_lengthen(L,M+2); }
    1019     2441182 :     L[m+1] = n;
    1020     2441182 :     if (n == 0) { if (++flag > 2 && m > Mestimate) break; } else flag = 0;
    1021             :   }
    1022       27820 :   m -= 2; while (m > 0 && !L[m]) m--;
    1023       26980 :   if (m == 0) { nmax = 1; L[1] = 1; m = 1; } /* can happen for tiny bitprec */
    1024       26980 :   setlg(L, m+1); S->M = m-1;
    1025       26980 :   S->L = L;
    1026       26980 :   S->nmax = nmax;
    1027             : 
    1028       26980 :   S->Dmax = S->D + (long)ceil((S->M * hd * S->MAXs - S->sub) / M_LN2);
    1029       26980 :   if (S->Dmax < S->D) S->Dmax = S->D;
    1030       26980 :   S->precmax = nbits2prec(S->Dmax);
    1031       26980 :   if (DEBUGLEVEL > 1)
    1032           0 :     err_printf("Dmax=%ld, D=%ld, M = %ld, nmax = %ld, m0 = %ld\n",
    1033             :                S->Dmax,S->D,S->M,S->nmax, S->m0);
    1034       26980 : }
    1035             : 
    1036             : static GEN
    1037       10878 : lfuninit_pol(GEN v, GEN poqk, long prec)
    1038             : {
    1039       10878 :   long m, M = lg(v) - 2;
    1040       10878 :   GEN pol = cgetg(M+3, t_POL);
    1041       10878 :   pol[1] = evalsigne(1) | evalvarn(0);
    1042       10878 :   gel(pol, 2) = gprec_w(gmul2n(gel(v,1), -1), prec);
    1043       10878 :   if (poqk)
    1044      523855 :     for (m = 2; m <= M+1; m++)
    1045      513033 :       gel(pol, m+1) = gprec_w(gmul(gel(poqk,m), gel(v,m)), prec);
    1046             :   else
    1047        2324 :     for (m = 2; m <= M+1; m++)
    1048        2268 :       gel(pol, m+1) = gprec_w(gel(v,m), prec);
    1049       10878 :   return RgX_renormalize_lg(pol, M+3);
    1050             : }
    1051             : 
    1052             : static void
    1053       78757 : worker_init(long q, GEN *an, GEN *bn, GEN *AB, GEN *A, GEN *B)
    1054             : {
    1055       78757 :   if (typ(*bn) == t_INT) *bn = NULL;
    1056       78757 :   if (*bn)
    1057             :   {
    1058         712 :     *AB = cgetg(3, t_VEC);
    1059         712 :     gel(*AB,1) = *A = cgetg(q+1, t_VEC);
    1060         712 :     gel(*AB,2) = *B = cgetg(q+1, t_VEC);
    1061         712 :     if (typ(an) == t_VEC) *an = RgV_kill0(*an);
    1062         712 :     if (typ(bn) == t_VEC) *bn = RgV_kill0(*bn);
    1063             :   }
    1064             :   else
    1065             :   {
    1066       78045 :     *B = NULL;
    1067       78045 :     *AB = *A = cgetg(q+1, t_VEC);
    1068       78049 :     if (typ(*an) == t_VEC) *an = RgV_kill0(*an);
    1069             :   }
    1070       78763 : }
    1071             : GEN
    1072       22393 : lfuninit_theta2_worker(long r, GEN L, GEN qk, GEN a, GEN di, GEN an, GEN bn)
    1073             : {
    1074       22393 :   long q, m, prec = di[1], M = di[2], m0 = di[3], L0 = lg(an)-1;
    1075             :   GEN AB, A, B;
    1076       22393 :   worker_init((M - r) / m0 + 1, &an, &bn, &AB, &A, &B);
    1077      302279 :   for (q = 0, m = r; m <= M; m += m0, q++)
    1078             :   {
    1079      279883 :     GEN t = gel(qk, m+1);
    1080      279883 :     long N = minss(L[m+1],L0);
    1081      279882 :     gel(A, q+1) = theta2(an, N, t, a, prec); /* theta(exp(mh)) */
    1082      279884 :     if (bn) gel(B, q+1) = theta2(bn, N, t, a, prec);
    1083             :   }
    1084       22396 :   return AB;
    1085             : }
    1086             : 
    1087             : /* theta(exp(mh)) ~ sum_{n <= N} a(n) k[m,n] */
    1088             : static GEN
    1089      239262 : an_msum(GEN an, long N, GEN vKm)
    1090             : {
    1091      239262 :   pari_sp av = avma;
    1092      239262 :   GEN s = gen_0;
    1093             :   long n;
    1094    12571430 :   for (n = 1; n <= N; n++)
    1095    12332744 :     if (gel(vKm,n))
    1096             :     {
    1097     6133827 :       GEN c = mul_an(an, n, gel(vKm,n));
    1098     6132901 :       if (c) s = gadd(s, c);
    1099             :     }
    1100      238686 :   return gc_upto(av, s);
    1101             : }
    1102             : 
    1103             : GEN
    1104       56380 : lfuninit_worker(long r, GEN K, GEN L, GEN peh2d, GEN vroots, GEN dr, GEN di,
    1105             :                 GEN an, GEN bn)
    1106             : {
    1107       56380 :   pari_sp av0 = avma;
    1108       56380 :   long m, n, q, L0 = lg(an)-1;
    1109       56380 :   double sig0 = rtodbl(gel(dr,1)), sub2 = rtodbl(gel(dr,2));
    1110       56373 :   double k1 = rtodbl(gel(dr,3)), MAXs = rtodbl(gel(dr,4));
    1111       56369 :   long D = di[1], M = di[2], m0 = di[3];
    1112       56369 :   double M0 = sig0? sub2 / sig0: 1./0.;
    1113       56369 :   GEN AB, A, B, vK = cgetg(M/m0 + 2, t_VEC);
    1114             : 
    1115      294581 :   for (q = 0, m = r; m <= M; m += m0, q++)
    1116      238221 :     gel(vK, q+1) = const_vec(L[m+1], NULL);
    1117       56360 :   worker_init(q, &an, &bn, &AB, &A, &B);
    1118      293785 :   for (m -= m0, q--; m >= 0; m -= m0, q--)
    1119             :   {
    1120      238271 :     double c1 = D + ((m > M0)? m * sig0 - sub2 : 0);
    1121      238271 :     GEN vKm = gel(vK,q+1); /* conceptually K(m,n) */
    1122    12578998 :     for (n = 1; n <= L[m+1]; n++)
    1123             :     {
    1124             :       GEN t2d, kmn;
    1125    12341580 :       long nn, mm, qq, p = 0;
    1126             :       double c, c2;
    1127             :       pari_sp av;
    1128             : 
    1129    12341580 :       if (gel(vKm, n)) continue; /* done already */
    1130     9259881 :       c = c1 + k1 * log2(n);
    1131             :       /* n *= 2; m -= m0 => c += c2, provided m >= M0. Else c += k1 */
    1132     9259881 :       c2 = k1 - MAXs;
    1133             :       /* p = largest (absolute) accuracy to which we need K(m,n) */
    1134    14692900 :       for (mm=m,nn=n; mm >= M0;)
    1135             :       {
    1136    11065747 :         if (nn <= L[mm+1] && (gel(an, nn) || (bn && gel(bn, nn))))
    1137     3002061 :           if (c > 0) p = maxuu(p, (ulong)c);
    1138    11066032 :         nn <<= 1;
    1139    11066032 :         mm -= m0; if (mm >= M0) c += c2; else { c += k1; break; }
    1140             :       }
    1141             :       /* mm < M0 || nn > L[mm+1] */
    1142    16510315 :       for (         ; mm >= 0; nn<<=1,mm-=m0,c+=k1)
    1143     7250240 :         if (nn <= L[mm+1] && (gel(an, nn) || (bn && gel(bn, nn))))
    1144     1771351 :           if (c > 0) p = maxuu(p, (ulong)c);
    1145     9260075 :       if (!p) continue; /* a_{n 2^v} = 0 for all v in range */
    1146     3057383 :       av = avma;
    1147     3057383 :       t2d = mpmul(gel(vroots,n), gel(peh2d,m+1));/*(n exp(mh)/sqrt(N))^(2/d)*/
    1148     3058090 :       kmn = gc_upto(av, gammamellininvrt(K, t2d, p));
    1149     9345169 :       for (qq=q,mm=m,nn=n; mm >= 0; nn<<=1,mm-=m0,qq--)
    1150     6288833 :         if (nn <= L[mm+1]) gmael(vK, qq+1, nn) = kmn;
    1151             :     }
    1152             :   }
    1153      293753 :   for (q = 0, m = r; m <= M; m += m0, q++)
    1154             :   {
    1155      238248 :     long N = minss(L0, L[m+1]);
    1156      238247 :     gel(A, q+1) = an_msum(an, N, gel(vK,q+1));
    1157      238239 :     if (bn) gel(B, q+1) = an_msum(bn, N, gel(vK,q+1));
    1158             :   }
    1159       55505 :   return gc_upto(av0, AB);
    1160             : }
    1161             : /* return A = [\theta(exp(mh)), m=0..M], theta(t) = sum a(n) K(n/sqrt(N) t),
    1162             :  * h = log(2)/m0. If bn != NULL, return the pair [A, B] */
    1163             : static GEN
    1164       10724 : lfuninit_ab(GEN theta, GEN h, struct lfunp *S)
    1165             : {
    1166       10724 :   const long M = S->M, prec = S->precmax;
    1167       10724 :   GEN tech = linit_get_tech(theta), isqN = theta_get_isqrtN(tech);
    1168       10724 :   GEN an = S->an, bn = S->bn, va, vb;
    1169             :   struct pari_mt pt;
    1170             :   GEN worker;
    1171             :   long m0, r, pending;
    1172             : 
    1173       10724 :   if (S->vgaell)
    1174             :   { /* d=2 and Vga = [a,a+1] */
    1175        7126 :     GEN a = vecmin(ldata_get_gammavec(linit_get_ldata(theta)));
    1176        7126 :     GEN qk = gpowers0(mpexp(h), M, isqN);
    1177        7126 :     m0 = minss(M+1, mt_nbthreads());
    1178        7126 :     worker = snm_closure(is_entry("_lfuninit_theta2_worker"),
    1179             :                          mkvecn(6, S->L, qk, a, mkvecsmall3(prec, M, m0),
    1180             :                                 an, bn? bn: gen_0));
    1181             :   }
    1182             :   else
    1183             :   {
    1184             :     GEN vroots, peh2d, d2;
    1185        3598 :     double sig0 = S->MAXs / S->m0, sub2 = S->sub / M_LN2;
    1186             :     /* For all 0<= m <= M, and all n <= L[m+1] such that a_n!=0, we compute
    1187             :      *   k[m,n] = K(n exp(mh)/sqrt(N))
    1188             :      * with ln(absolute error) <= E + max(mh sigma - sub, 0) + k1 * log(n).
    1189             :      * N.B. we use the 'rt' variant and pass (n exp(mh)/sqrt(N))^(2/d).
    1190             :      * Speedup: if n' = 2n and m' = m - m0 >= 0; then k[m,n] = k[m',n']. */
    1191        3598 :     vroots = mkvroots(S->d, S->nmax, prec); /* vroots[n] = n^(2/d) */
    1192        3598 :     d2 = gdivgu(gen_2, S->d);
    1193        3598 :     peh2d = gpowers0(gexp(gmul(d2,h), prec), M, gpow(isqN, d2, prec));
    1194        3598 :     m0 = S->m0; /* peh2d[m+1] = (exp(mh)/sqrt(N))^(2/d) */
    1195        3598 :     worker = snm_closure(is_entry("_lfuninit_worker"),
    1196             :                          mkvecn(8, theta_get_K(tech), S->L, peh2d, vroots,
    1197             :                                 mkvec4(dbltor(sig0), dbltor(sub2),
    1198             :                                        dbltor(S->k1), dbltor(S->MAXs)),
    1199             :                                 mkvecsmall3(S->D, M, m0),
    1200             :                                 an, bn? bn: gen_0));
    1201             :     /* For each 0 <= m <= M, we will sum for n<=L[m+1] a(n) K(m,n)
    1202             :      * bit accuracy for K(m,n): D + k1*log2(n) + 1_{m > M0} (m*sig0 - sub2)
    1203             :      * We restrict m to arithmetic progressions r mod m0 to save memory and
    1204             :      * allow parallelization */
    1205             :   }
    1206       10724 :   va = cgetg(M+2, t_VEC);
    1207       10724 :   vb = bn? cgetg(M+2, t_VEC): NULL;
    1208       10724 :   mt_queue_start_lim(&pt, worker, m0);
    1209       10724 :   pending = 0;
    1210      110733 :   for (r = 0; r < m0 || pending; r++)
    1211             :   { /* m = q m0 + r */
    1212             :     GEN done, A, B;
    1213             :     long q, m, workid;
    1214      100009 :     mt_queue_submit(&pt, r, r < m0 ? mkvec(utoi(r)): NULL);
    1215      100009 :     done = mt_queue_get(&pt, &workid, &pending);
    1216      100009 :     if (!done) continue;
    1217       78782 :     if (bn) { A = gel(done,1); B = gel(done,2); } else { A = done; B = NULL; }
    1218      596981 :     for (q = 0, m = workid; m <= M; m += m0, q++)
    1219             :     {
    1220      518199 :       gel(va, m+1) = gel(A, q+1);
    1221      518199 :       if (bn) gel(vb, m+1) = gel(B, q+1);
    1222             :     }
    1223             :   }
    1224       10724 :   mt_queue_end(&pt);
    1225       10724 :   return bn? mkvec2(va, vb): va;
    1226             : }
    1227             : 
    1228             : static void
    1229      140832 : parse_dom(double k, GEN dom, struct lfunp *S)
    1230             : {
    1231      140832 :   long l = lg(dom);
    1232      140832 :   if (typ(dom)!=t_VEC) pari_err_TYPE("lfuninit [domain]", dom);
    1233      140832 :   if (l == 1)
    1234             :   {
    1235          98 :     S->dc = 0;
    1236          98 :     S->dw = -1;
    1237          98 :     S->dh = -1; return;
    1238             :   }
    1239      140734 :   if (l == 2)
    1240             :   {
    1241       38124 :     S->dc = k/2.;
    1242       38124 :     S->dw = 0.;
    1243       38124 :     S->dh = gtodouble(gel(dom,1));
    1244             :   }
    1245      102610 :   else if (l == 3)
    1246             :   {
    1247         301 :     S->dc = k/2.;
    1248         301 :     S->dw = gtodouble(gel(dom,1));
    1249         301 :     S->dh = gtodouble(gel(dom,2));
    1250             :   }
    1251      102309 :   else if (l == 4)
    1252             :   {
    1253      102309 :     S->dc = gtodouble(gel(dom,1));
    1254      102309 :     S->dw = gtodouble(gel(dom,2));
    1255      102309 :     S->dh = gtodouble(gel(dom,3));
    1256             :   }
    1257             :   else
    1258             :   {
    1259           0 :     pari_err_TYPE("lfuninit [domain]", dom);
    1260           0 :     S->dc = S->dw = S->dh = 0; /*-Wall*/
    1261             :   }
    1262      140734 :   if (S->dw < 0 || S->dh < 0) pari_err_TYPE("lfuninit [domain]", dom);
    1263             : }
    1264             : 
    1265             : /* do we have dom \subset dom0 ? dom = [center, width, height] */
    1266             : int
    1267       25029 : sdomain_isincl(double k, GEN dom, GEN dom0)
    1268             : {
    1269             :   struct lfunp S0, S;
    1270       25029 :   parse_dom(k, dom, &S); if (S.dw < 0) return 1;
    1271       25029 :   parse_dom(k, dom0, &S0); if (S0.dw < 0) return 0;
    1272       25029 :   return S0.dc - S0.dw <= S.dc - S.dw
    1273       25029 :       && S0.dc + S0.dw >= S.dc + S.dw && S0.dh >= S.dh;
    1274             : }
    1275             : 
    1276             : static int
    1277       25106 : checklfuninit(GEN linit, GEN DOM, long der, long bitprec)
    1278             : {
    1279       25106 :   GEN ldata = linit_get_ldata(linit);
    1280       25106 :   GEN domain = lfun_get_domain(linit_get_tech(linit));
    1281       25106 :   GEN dom = domain_get_dom(domain);
    1282       25106 :   if (lg(dom) == 1) return 1;
    1283       25029 :   return domain_get_der(domain) >= der
    1284       25029 :     && domain_get_bitprec(domain) >= bitprec
    1285       50058 :     && sdomain_isincl(gtodouble(ldata_get_k(ldata)), DOM, dom);
    1286             : }
    1287             : 
    1288             : static GEN
    1289        2387 : ginvsqrtvec(GEN x, long prec)
    1290             : {
    1291        2387 :   if (is_vec_t(typ(x)))
    1292        1813 :     pari_APPLY_same(ginv(gsqrt(gel(x,i), prec)))
    1293        1939 :   else return ginv(gsqrt(x, prec));
    1294             : }
    1295             : 
    1296             : GEN
    1297       11830 : lfuninit_make(long t, GEN ldata, GEN tech, GEN domain)
    1298             : {
    1299       11830 :   GEN Vga = ldata_get_gammavec(ldata);
    1300       11830 :   long d = lg(Vga)-1;
    1301       11830 :   GEN w2 = gen_1, k2 = gmul2n(ldata_get_k(ldata), -1);
    1302       11830 :   GEN expot = gdivgu(gadd(gmulsg(d, gsubgs(k2, 1)), sumVga(Vga)), 4);
    1303       11830 :   if (typ(ldata_get_dual(ldata))==t_INT)
    1304             :   {
    1305       11676 :     GEN eno = ldata_get_rootno(ldata);
    1306       11676 :     long prec = nbits2prec( domain_get_bitprec(domain) );
    1307       11676 :     if (!isint1(eno)) w2 = ginvsqrtvec(eno, prec);
    1308             :   }
    1309       11830 :   tech = mkvec3(domain, tech, mkvec4(k2, w2, expot, gammafactor(Vga)));
    1310       11830 :   return mkvec3(mkvecsmall(t), ldata, tech);
    1311             : }
    1312             : static GEN
    1313         224 : lfunnoinit(GEN ldata, long bitprec)
    1314             : {
    1315         224 :   GEN tech, domain = mkvec2(cgetg(1,t_VEC), mkvecsmall2(0, bitprec));
    1316         224 :   GEN R = gen_0, r = ldata_get_residue(ldata), v = lfunrootres(ldata, bitprec);
    1317         224 :   ldata = shallowcopy(ldata);
    1318         224 :   gel(ldata,6) = gel(v,3);
    1319         224 :   if (r)
    1320             :   {
    1321         196 :     if (isintzero(r)) setlg(ldata,7); else gel(ldata,7) = r;
    1322         196 :     R = gel(v,2);
    1323             :   }
    1324         224 :   tech = mkvec3(domain, gen_0, R);
    1325         224 :   return lfuninit_make(t_LDESC_INIT, ldata, tech, domain);
    1326             : }
    1327             : 
    1328             : static void
    1329        3598 : lfunparams2(struct lfunp *S)
    1330             : {
    1331        3598 :   GEN L = S->L, an = S->an, bn = S->bn;
    1332             :   double pmax;
    1333        3598 :   long m, nan, nmax, neval, M = S->M;
    1334             : 
    1335        3598 :   S->vgaell = 0;
    1336             :   /* try to reduce parameters now we know the a_n (some may be 0) */
    1337        3598 :   if (typ(an) == t_VEC) an = RgV_kill0(an);
    1338        3598 :   nan = S->nmax; /* lg(an)-1 may be large than this */
    1339        3598 :   nmax = neval = 0;
    1340        3598 :   if (!bn)
    1341      240866 :     for (m = 0; m <= M; m++)
    1342             :     {
    1343      237289 :       long n = minss(nan, L[m+1]);
    1344      341897 :       while (n > 0 && !gel(an,n)) n--;
    1345      237289 :       if (n > nmax) nmax = n;
    1346      237289 :       neval += n;
    1347      237289 :       L[m+1] = n; /* reduce S->L[m+1] */
    1348             :     }
    1349             :   else
    1350             :   {
    1351          21 :     if (typ(bn) == t_VEC) bn = RgV_kill0(bn);
    1352        1036 :     for (m = 0; m <= M; m++)
    1353             :     {
    1354        1015 :       long n = minss(nan, L[m+1]);
    1355        1057 :       while (n > 0 && !gel(an,n) && !gel(bn,n)) n--;
    1356        1015 :       if (n > nmax) nmax = n;
    1357        1015 :       neval += n;
    1358        1015 :       L[m+1] = n; /* reduce S->L[m+1] */
    1359             :     }
    1360             :   }
    1361        3598 :   if (DEBUGLEVEL >= 1) err_printf("expected evaluations: %ld\n", neval);
    1362        3598 :   for (; M > 0; M--)
    1363        3598 :     if (L[M+1]) break;
    1364        3598 :   setlg(L, M+2);
    1365        3598 :   S->M = M;
    1366        3598 :   S->nmax = nmax;
    1367             : 
    1368             :   /* need K(n*exp(mh)/sqrt(N)) to absolute accuracy
    1369             :    *   D + k1*log(n) + max(m * sig0 - sub2, 0) */
    1370        3598 :   pmax = S->D + S->k1 * log2(L[1]);
    1371        3598 :   if (S->MAXs)
    1372             :   {
    1373        3598 :     double sig0 = S->MAXs/S->m0, sub2 = S->sub / M_LN2;
    1374      202484 :     for (m = ceil(sub2 / sig0); m <= S->M; m++)
    1375             :     {
    1376      198886 :       double c = S->D + m*sig0 - sub2;
    1377      198886 :       if (S->k1 > 0) c += S->k1 * log2(L[m+1]);
    1378      198886 :       pmax = maxdd(pmax, c);
    1379             :     }
    1380             :   }
    1381        3598 :   S->Dmax = pmax;
    1382        3598 :   S->precmax = nbits2prec(pmax);
    1383        3598 : }
    1384             : 
    1385             : static GEN
    1386       10738 : lfun_init_theta(GEN ldata, GEN eno, struct lfunp *S)
    1387             : {
    1388       10738 :   GEN an2, dual, tdom = NULL, Vga = ldata_get_gammavec(ldata);
    1389       10738 :   long L, extrabit = 0, prec = S->precmax;
    1390       10738 :   if (eno)
    1391        6524 :     L = S->nmax;
    1392             :   else
    1393             :   {
    1394        4214 :     tdom = dbltor(sqrt(0.5));
    1395        4214 :     L = maxss(S->nmax, lfunthetacost(ldata, tdom, 0, S->D, &extrabit));
    1396        4214 :     prec += nbits2extraprec(extrabit);
    1397             :   }
    1398       10738 :   dual = ldata_get_dual(ldata);
    1399       10738 :   S->an = ldata_vecan(ldata_get_an(ldata), L, prec);
    1400       10724 :   S->bn = typ(dual)==t_INT? NULL: ldata_vecan(dual, S->nmax, prec);
    1401       10724 :   if (!vgaell(Vga)) lfunparams2(S);
    1402             :   else
    1403             :   {
    1404        7126 :     S->an = antwist(S->an, Vga, prec);
    1405        7126 :     if (S->bn) S->bn = antwist(S->bn, Vga, prec);
    1406        7126 :     S->vgaell = 1;
    1407             :   }
    1408       10724 :   an2 = lg(Vga)-1 == 1? antwist(S->an, Vga, prec): S->an;
    1409       10724 :   return lfunthetainit0(ldata, tdom, an2, 0, S->Dmax, extrabit);
    1410             : }
    1411             : 
    1412             : GEN
    1413       16242 : lfuncost(GEN L, GEN dom, long der, long bit)
    1414             : {
    1415       16242 :   pari_sp av = avma;
    1416       16242 :   GEN ldata = lfunmisc_to_ldata_shallow(L);
    1417       16242 :   GEN w, k = ldata_get_k(ldata);
    1418             :   struct lfunp S;
    1419             : 
    1420       16242 :   parse_dom(gtodouble(k), dom, &S); if (S.dw < 0) return mkvecsmall2(0, 0);
    1421       16242 :   lfunp_set(ldata, der, bit, &S);
    1422       16242 :   w = ldata_get_rootno(ldata);
    1423       16242 :   if (isintzero(w)) /* for lfunrootres */
    1424           7 :     S.nmax = maxss(S.nmax, lfunthetacost(ldata,dbltor(sqrt(0.5)),0,bit+1,NULL));
    1425       16242 :   set_avma(av); return mkvecsmall2(S.nmax, S.Dmax);
    1426             : }
    1427             : GEN
    1428          49 : lfuncost0(GEN L, GEN dom, long der, long bitprec)
    1429             : {
    1430          49 :   pari_sp av = avma;
    1431             :   GEN C;
    1432             : 
    1433          49 :   if (is_linit(L))
    1434             :   {
    1435          28 :     GEN tech = linit_get_tech(L);
    1436          28 :     GEN domain = lfun_get_domain(tech);
    1437          28 :     dom = domain_get_dom(domain);
    1438          28 :     der = domain_get_der(domain);
    1439          28 :     bitprec = domain_get_bitprec(domain);
    1440          28 :     if (linit_get_type(L) == t_LDESC_PRODUCT)
    1441             :     {
    1442          21 :       GEN v = lfunprod_get_fact(linit_get_tech(L)), F = gel(v,1);
    1443          21 :       long i, l = lg(F);
    1444          21 :       C = cgetg(l, t_VEC);
    1445          70 :       for (i = 1; i < l; ++i)
    1446          49 :         gel(C, i) = zv_to_ZV( lfuncost(gel(F,i), dom, der, bitprec) );
    1447          21 :       return gc_upto(av, C);
    1448             :     }
    1449             :   }
    1450          28 :   if (!dom) pari_err_TYPE("lfuncost [missing s domain]", L);
    1451          28 :   C = lfuncost(L,dom,der,bitprec);
    1452          28 :   return gc_upto(av, zv_to_ZV(C));
    1453             : }
    1454             : 
    1455             : static int
    1456        9716 : is_dirichlet(GEN ldata)
    1457             : {
    1458        9716 :   switch(ldata_get_type(ldata))
    1459             :   {
    1460        1330 :     case t_LFUN_ZETA:
    1461             :     case t_LFUN_KRONECKER:
    1462        1330 :     case t_LFUN_CHIZ: return 1;
    1463         980 :     case t_LFUN_CHIGEN: return ldata_get_degree(ldata)==1;
    1464        7406 :     default: return 0;
    1465             :   }
    1466             : }
    1467             : 
    1468             : static ulong
    1469       11016 : lfuninit_cutoff(GEN ldata)
    1470             : {
    1471       11016 :   GEN gN = ldata_get_conductor(ldata);
    1472             :   ulong L, N;
    1473       11016 :   if (ldata_get_type(ldata) == t_LFUN_NF) /* N ~ f^(d-1), exact for d prime */
    1474         742 :     gN = sqrtnint(gN, ldata_get_degree(ldata) - 1);
    1475       11016 :   N = itou_or_0(gN);
    1476       11016 :   if (N > 1000) L = 7000;
    1477       11002 :   else if (N > 100) L = 5000;
    1478        8013 :   else if (N > 15) L = 3000;
    1479        7558 :   else L = 2500;
    1480       11016 :   return L;
    1481             : }
    1482             : 
    1483             : GEN
    1484       36642 : lfuninit(GEN lmisc, GEN dom, long der, long bitprec)
    1485             : {
    1486       36642 :   pari_sp av = avma;
    1487             :   GEN poqk, AB, R, h, theta, ldata, eno, r, domain, tech, k;
    1488             :   struct lfunp S;
    1489             : 
    1490       36642 :   if (is_linit(lmisc))
    1491             :   {
    1492       25155 :     long t = linit_get_type(lmisc);
    1493       25155 :     if (t == t_LDESC_INIT || t == t_LDESC_PRODUCT)
    1494             :     {
    1495       25106 :       if (checklfuninit(lmisc, dom, der, bitprec)) return lmisc;
    1496          21 :       pari_warn(warner,"lfuninit: insufficient initialization");
    1497             :     }
    1498             :   }
    1499       11557 :   ldata = lfunmisc_to_ldata_shallow(lmisc);
    1500             : 
    1501       11557 :   switch (ldata_get_type(ldata))
    1502             :   {
    1503         630 :   case t_LFUN_NF:
    1504             :     {
    1505         630 :       GEN T = gel(ldata_get_an(ldata), 2);
    1506         630 :       return gc_GEN(av, lfunzetakinit(T, dom, der, bitprec));
    1507             :     }
    1508          91 :   case t_LFUN_ABELREL:
    1509             :     {
    1510          91 :       GEN T = gel(ldata_get_an(ldata), 2);
    1511          91 :       return gc_GEN(av, lfunabelianrelinit(gel(T,1), gel(T,2), dom, der, bitprec));
    1512             :     }
    1513             :   }
    1514       10836 :   k = ldata_get_k(ldata);
    1515       10836 :   parse_dom(gtodouble(k), dom, &S);
    1516             :   /* Reduce domain for Dirichlet characters. NOT for Abelian t_LFUN_NF,
    1517             :    * handled above. */
    1518       10836 :   if (S.dw >= 0 && (!der && is_dirichlet(ldata)))
    1519        1750 :     S.dh = mindd(S.dh, lfuninit_cutoff(ldata));
    1520       10836 :   if (S.dw < 0)
    1521             :   {
    1522          98 :     if (der)
    1523           0 :       pari_err_IMPL("domain = [] for derivatives in lfuninit");
    1524          98 :     if (!is_dirichlet(ldata))
    1525           0 :       pari_err_IMPL("domain = [] for L functions of degree > 1");
    1526          98 :     return gc_GEN(av, lfunnoinit(ldata, bitprec));
    1527             :   }
    1528             : 
    1529       10738 :   lfunp_set(ldata, der, bitprec, &S);
    1530       10738 :   ldata = ldata_newprec(ldata, nbits2prec(S.Dmax));
    1531       10738 :   r = ldata_get_residue(ldata);
    1532             :   /* Note: all guesses should already have been performed (thetainit more
    1533             :    * expensive than needed: should be either tdom = 1 or bitprec = S.D).
    1534             :    * BUT if the root number / polar part do not have an algebraic
    1535             :    * expression, there is no way to do this until we know the
    1536             :    * precision, i.e. now. So we can't remove guessing code from here and
    1537             :    * lfun_init_theta */
    1538       10738 :   if (r && isintzero(r)) eno = NULL;
    1539             :   else
    1540             :   {
    1541       10738 :     eno = ldata_get_rootno(ldata);
    1542       10738 :     if (isintzero(eno)) eno = NULL;
    1543             :   }
    1544       10738 :   theta = lfun_init_theta(ldata, eno, &S);
    1545       10724 :   if (eno && !r)
    1546        4599 :     R = gen_0;
    1547             :   else
    1548             :   {
    1549        6125 :     GEN v = lfunrootres(theta, S.D);
    1550        6125 :     ldata = shallowcopy(ldata);
    1551        6125 :     gel(ldata, 6) = gel(v,3);
    1552        6125 :     r = gel(v,1);
    1553        6125 :     R = gel(v,2);
    1554        6125 :     if (isintzero(r)) setlg(ldata,7); else gel(ldata, 7) = r;
    1555             :   }
    1556       10724 :   h = divru(mplog2(S.precmax), S.m0);
    1557             :   /* exp(kh/2 . [0..M]) */
    1558       10724 :   poqk = gequal0(k) ? NULL
    1559       10724 :        : gpowers(gprec_w(mpexp(gmul2n(gmul(k,h), -1)), S.precmax), S.M);
    1560       10724 :   AB = lfuninit_ab(theta, h, &S);
    1561       10724 :   if (S.bn)
    1562             :   {
    1563         154 :     GEN A = gel(AB,1), B = gel(AB,2);
    1564         154 :     A = lfuninit_pol(A, poqk, S.precmax);
    1565         154 :     B = lfuninit_pol(B, poqk, S.precmax);
    1566         154 :     AB = mkvec2(A, B);
    1567             :   }
    1568             :   else
    1569       10570 :     AB = lfuninit_pol(AB, poqk, S.precmax);
    1570       10724 :   tech = mkvec3(h, AB, R);
    1571       10724 :   domain = mkvec2(dom, mkvecsmall2(der, bitprec));
    1572       10724 :   return gc_GEN(av, lfuninit_make(t_LDESC_INIT, ldata, tech, domain));
    1573             : }
    1574             : 
    1575             : GEN
    1576         532 : lfuninit0(GEN lmisc, GEN dom, long der, long bitprec)
    1577             : {
    1578         532 :   GEN z = lfuninit(lmisc, dom, der, bitprec);
    1579         532 :   return z == lmisc? gcopy(z): z;
    1580             : }
    1581             : 
    1582             : /* If s is a pole of Lambda, return polar part at s; else return NULL */
    1583             : static GEN
    1584        5151 : lfunpoleresidue(GEN R, GEN s)
    1585             : {
    1586             :   long j;
    1587       14676 :   for (j = 1; j < lg(R); j++)
    1588             :   {
    1589       10085 :     GEN Rj = gel(R, j), be = gel(Rj, 1);
    1590       10085 :     if (gequal(s, be)) return gel(Rj, 2);
    1591             :   }
    1592        4591 :   return NULL;
    1593             : }
    1594             : 
    1595             : /* Compute contribution of polar part at s when not a pole. */
    1596             : static GEN
    1597        8377 : veccothderivn(GEN a, long n)
    1598             : {
    1599             :   long i;
    1600        8377 :   pari_sp av = avma;
    1601        8377 :   GEN c = pol_x(0), cp = mkpoln(3, gen_m1, gen_0, gen_1);
    1602        8377 :   GEN v = cgetg(n+2, t_VEC);
    1603        8377 :   gel(v, 1) = poleval(c, a);
    1604       25250 :   for(i = 2; i <= n+1; i++)
    1605             :   {
    1606       16873 :     c = ZX_mul(ZX_deriv(c), cp);
    1607       16873 :     gel(v, i) = gdiv(poleval(c, a), mpfact(i-1));
    1608             :   }
    1609        8377 :   return gc_GEN(av, v);
    1610             : }
    1611             : 
    1612             : static GEN
    1613        8496 : polepart(long n, GEN h, GEN C)
    1614             : {
    1615        8496 :   GEN h2n = gpowgs(gdiv(h, gen_2), n-1);
    1616        8496 :   GEN res = gmul(h2n, gel(C,n));
    1617        8496 :   return odd(n)? res : gneg(res);
    1618             : }
    1619             : 
    1620             : static GEN
    1621        4157 : lfunsumcoth(GEN R, GEN s, GEN h, long prec)
    1622             : {
    1623             :   long i,j;
    1624        4157 :   GEN S = gen_0;
    1625       12534 :   for (j = 1; j < lg(R); ++j)
    1626             :   {
    1627        8377 :     GEN r = gel(R,j), be = gel(r,1), Rj = gel(r, 2);
    1628        8377 :     long e = valser(Rj);
    1629        8377 :     GEN z1 = gexpm1(gmul(h, gsub(s,be)), prec); /* exp(h(s-beta))-1 */
    1630        8377 :     GEN c1 = gaddgs(gdivsg(2, z1), 1); /* coth((h/2)(s-beta)) */
    1631        8377 :     GEN C1 = veccothderivn(c1, 1-e);
    1632       16873 :     for (i = e; i < 0; i++)
    1633             :     {
    1634        8496 :       GEN Rbe = mysercoeff(Rj, i);
    1635        8496 :       GEN p1 = polepart(-i, h, C1);
    1636        8496 :       S = gadd(S, gmul(Rbe, p1));
    1637             :     }
    1638             :   }
    1639        4157 :   return gmul2n(S, -1);
    1640             : }
    1641             : 
    1642             : static GEN lfunlambda_OK(GEN linit, GEN s, GEN sdom, long bitprec);
    1643             : /* L is a t_LDESC_PRODUCT or t_LDESC_INIT Linit */
    1644             : static GEN
    1645        2456 : _product(GEN (*fun)(GEN,GEN,long), GEN L, GEN s, long bitprec)
    1646             : {
    1647        2456 :   GEN ldata = linit_get_ldata(L), v, r, F, E, C, cs;
    1648             :   long i, l;
    1649             :   int isreal;
    1650        2456 :   if (linit_get_type(L) == t_LDESC_INIT) return fun(ldata, s, bitprec);
    1651        1924 :   v = lfunprod_get_fact(linit_get_tech(L));
    1652        1924 :   F = gel(v,1); E = gel(v,2); C = gel(v,3); l = lg(F);
    1653        1924 :   cs = conj_i(s); isreal = gequal(imag_i(s), imag_i(cs));
    1654        6374 :   for (i = 1, r = gen_1; i < l; ++i)
    1655             :   {
    1656        4450 :     GEN f = fun(gel(F, i), s, bitprec);
    1657        4450 :     if (typ(f)==t_VEC) f = RgV_prod(f);
    1658        4450 :     if (E[i]) r = gmul(r, gpowgs(f, E[i]));
    1659        4450 :     if (C[i])
    1660             :     {
    1661           0 :       GEN fc = isreal? f: conj_i(fun(gel(F, i), cs, bitprec));
    1662           0 :       r = gmul(r, gpowgs(fc, C[i]));
    1663             :     }
    1664             :   }
    1665        1924 :   return (ldata_isreal(ldata) && gequal0(imag_i(s)))? real_i(r): r;
    1666             : }
    1667             : 
    1668             : /* s a t_SER; # terms - 1 */
    1669             : static long
    1670        2248 : der_level(GEN s)
    1671        2248 : { return signe(s)? lg(s)-3: valser(s)-1; }
    1672             : 
    1673             : /* s a t_SER; return coeff(s, X^0) */
    1674             : static GEN
    1675        1232 : ser_coeff0(GEN s) { return simplify_shallow(polcoef_i(s, 0, -1)); }
    1676             : 
    1677             : static GEN
    1678       19850 : get_domain(GEN s, GEN *dom, long *der)
    1679             : {
    1680       19850 :   GEN sa = s;
    1681       19850 :   *der = 0;
    1682       19850 :   switch(typ(s))
    1683             :   {
    1684           7 :     case t_POL:
    1685           7 :     case t_RFRAC: s = toser_i(s);
    1686        1232 :     case t_SER:
    1687        1232 :       *der = der_level(s);
    1688        1232 :       sa = ser_coeff0(s);
    1689             :   }
    1690       19850 :   *dom = mkvec3(real_i(sa), gen_0, gabs(imag_i(sa),DEFAULTPREC));
    1691       19850 :   return s;
    1692             : }
    1693             : /* assume s went through get_domain and s/bitprec belong to domain */
    1694             : static GEN
    1695       33527 : lfunlambda_OK(GEN linit, GEN s, GEN sdom, long bitprec)
    1696             : {
    1697       33527 :   GEN dom, eno, ldata, tech, h, pol, k2, cost, S, S0 = NULL;
    1698             :   long prec, prec0;
    1699             :   struct lfunp D, D0;
    1700             : 
    1701       33527 :   if (linit_get_type(linit) == t_LDESC_PRODUCT)
    1702        1679 :     return _product(&lfunlambda, linit, s, bitprec);
    1703       31848 :   ldata = linit_get_ldata(linit);
    1704       31848 :   eno = ldata_get_rootno(ldata);
    1705       31848 :   tech = linit_get_tech(linit);
    1706       31848 :   dom = lfun_get_dom(tech);
    1707       31848 :   if (lg(dom) == 1) return lfunlambda(linit, s, bitprec); /* FIXME:not OK! */
    1708       31848 :   h = lfun_get_step(tech); prec = realprec(h);
    1709             :   /* try to reduce accuracy */
    1710       31848 :   parse_dom(0, sdom, &D0);
    1711       31848 :   parse_dom(0, dom, &D);
    1712       31848 :   if (0.8 * D.dh > D0.dh)
    1713             :   {
    1714       16165 :     cost = lfuncost(linit, sdom, typ(s)==t_SER? der_level(s): 0, bitprec);
    1715       16165 :     prec0 = nbits2prec(cost[2]);
    1716       16165 :     if (prec0 < prec) { prec = prec0; h = gprec_w(h, prec); }
    1717             :   }
    1718       31848 :   pol = lfun_get_pol(tech);
    1719       31848 :   s = gprec_w(s, prec);
    1720       31848 :   if (ldata_get_residue(ldata))
    1721             :   {
    1722        4556 :     GEN R = lfun_get_Residue(tech);
    1723        4556 :     GEN Ra = lfunpoleresidue(R, s);
    1724        4556 :     if (Ra) return gprec_w(Ra, nbits2prec(bitprec));
    1725        4157 :     S0 = lfunsumcoth(R, s, h, prec);
    1726             :   }
    1727       31449 :   k2 = lfun_get_k2(tech);
    1728       31449 :   if (typ(pol)==t_POL && typ(s) != t_SER && gequal(real_i(s), k2))
    1729       25030 :   { /* on critical line: shortcut */
    1730       25030 :     GEN polz, b = imag_i(s);
    1731       25030 :     polz = gequal0(b)? poleval(pol,gen_1): poleval(pol, expIr(gmul(h,b)));
    1732       25030 :     S = gadd(polz, gmulvec(eno, conj_i(polz)));
    1733             :   }
    1734             :   else
    1735             :   {
    1736        6419 :     GEN z = gexp(gmul(h, gsub(s, k2)), prec);
    1737        6419 :     GEN zi = ginv(z), zc = conj_i(zi);
    1738        6419 :     if (typ(pol)==t_POL)
    1739        6223 :       S = gadd(poleval(pol, z), gmulvec(eno, conj_i(poleval(pol, zc))));
    1740             :     else
    1741         196 :       S = gadd(poleval(gel(pol,1), z), gmulvec(eno, poleval(gel(pol,2), zi)));
    1742             :   }
    1743       31449 :   if (S0) S = gadd(S,S0);
    1744       31449 :   return gprec_w(gmul(S,h), nbits2prec(bitprec));
    1745             : }
    1746             : 
    1747             : static long
    1748       38308 : lfunspec_OK(GEN lmisc, GEN s, GEN *pldata)
    1749             : {
    1750       38308 :   long t, large = 0;
    1751             :   GEN ldata;
    1752       38308 :   *pldata = ldata = lfunmisc_to_ldata_shallow(lmisc);
    1753       38301 :   if (!is_linit(lmisc)) lmisc = ldata;
    1754       25841 :   else switch(linit_get_type(lmisc))
    1755             :   {
    1756       25799 :     case t_LDESC_INIT: case t_LDESC_PRODUCT:
    1757       25799 :       if (lg(lfun_get_dom(linit_get_tech(lmisc))) == 1) large = 1;
    1758       25799 :       break;
    1759          42 :     default: return 0;
    1760             :   }
    1761       38259 :   switch(typ(s))
    1762             :   {
    1763       37916 :     case t_INT: case t_REAL: case t_FRAC: case t_COMPLEX: break;
    1764         343 :     default: return 0;
    1765             :   }
    1766       37916 :   t = ldata_get_type(ldata);
    1767       37916 :   switch(t)
    1768             :   {
    1769        7930 :     case t_LFUN_KRONECKER: case t_LFUN_ZETA:
    1770        7930 :       if (typ(s) == t_INT && !is_bigint(s)) return 1;
    1771             :       /* fall through */
    1772             :     case t_LFUN_NF: case t_LFUN_CHIZ:
    1773        8056 :       if (!large)
    1774        7762 :         large = (fabs(gtodouble(imag_i(s))) >= lfuninit_cutoff(ldata));
    1775        8056 :       break;
    1776        2057 :     case t_LFUN_CHIGEN:
    1777        2057 :       if (ldata_get_degree(ldata) != 1) return 0;
    1778        1735 :       if (!large)
    1779        1483 :         large = (fabs(gtodouble(imag_i(s))) >= lfuninit_cutoff(ldata));
    1780        1735 :       break;
    1781             :   }
    1782       33513 :   if (large)
    1783             :   {
    1784        1904 :     if (t == t_LFUN_NF)
    1785             :     {
    1786          14 :       GEN an = ldata_get_an(ldata), nf = gel(an,2), G = galoisinit(nf, NULL);
    1787          14 :       if (isintzero(G) || !group_isabelian(galois_group(G))) return 0;
    1788             :     }
    1789        1904 :     return 2;
    1790             :   }
    1791       31609 :   return 0;
    1792             : }
    1793             : 
    1794             : GEN
    1795        5094 : lfunlambda(GEN lmisc, GEN s, long bitprec)
    1796             : {
    1797        5094 :   pari_sp av = avma;
    1798        5094 :   GEN linit = NULL, dom, z;
    1799             :   long der;
    1800        5094 :   s = get_domain(s, &dom, &der);
    1801        5094 :   if (!der)
    1802             :   {
    1803             :     GEN ldata;
    1804        4247 :     long t = lfunspec_OK(lmisc, s, &ldata);
    1805        4247 :     if (t == 1)
    1806             :     { /* special value ? */
    1807         553 :       GEN z = lfun(ldata, s, bitprec), gv = ldata_get_gammavec(ldata);
    1808         553 :       long e = itou(gel(gv, 1));
    1809         553 :       if (!isintzero(z) && (e || gsigne(s) > 0)) /* TODO */
    1810             :       {
    1811         476 :         GEN q = ldata_get_conductor(ldata);
    1812         476 :         long prec = nbits2prec(bitprec);
    1813         476 :         GEN se, r, Q = divir(q, mppi(prec));
    1814         476 :         se = gmul2n(gaddgs(s, e), -1);
    1815         476 :         r = gmul(gpow(Q, se, prec), ggamma(se, prec));
    1816         476 :         if (e && !equali1(q)) r = gdiv(r, gsqrt(q, prec));
    1817         504 :         return gc_upto(av, gmul(r, z));
    1818             :       }
    1819             :     }
    1820        3771 :     if (is_linit(lmisc)) linit = lmisc; else lmisc = ldata;
    1821        3771 :     if (t == 2)
    1822          28 :       return gc_GEN(av, linit? _product(&lfunlambda, linit, s, bitprec)
    1823           0 :                                    : lfunlambdalarge(ldata, s, bitprec));
    1824             :   }
    1825        4590 :   linit = lfuninit(lmisc, dom, der, bitprec);
    1826        4590 :   z = lfunlambda_OK(linit,s, dom, bitprec);
    1827        4590 :   return gc_GEN(av, z);
    1828             : }
    1829             : 
    1830             : static long
    1831       18459 : is_ser(GEN x)
    1832             : {
    1833       18459 :   long t = typ(x);
    1834       18459 :   if (t == t_SER) return 1;
    1835       16478 :   if (!is_vec_t(t) || lg(x)==1) return 0;
    1836         350 :   if (typ(gel(x,1))==t_SER) return 1;
    1837         252 :   return 0;
    1838             : }
    1839             : 
    1840             : static GEN
    1841         371 : lfunser(GEN L)
    1842             : {
    1843         371 :   long v = valser(L);
    1844         371 :   if (v > 0) return gen_0;
    1845         329 :   if (v == 0) L = gel(L, 2);
    1846             :   else
    1847         203 :     setlg(L, minss(lg(L), 2-v));
    1848         329 :   return L;
    1849             : }
    1850             : 
    1851             : static GEN
    1852         371 : lfunservec(GEN x)
    1853             : {
    1854         371 :   if (typ(x)==t_SER) return lfunser(x);
    1855           0 :   pari_APPLY_same(lfunser(gel(x,i)))
    1856             : }
    1857             : static GEN
    1858         105 : lfununext(GEN L)
    1859             : {
    1860         105 :   setlg(L, maxss(lg(L)-1, valser(L)? 2: 3));
    1861         105 :   return normalizeser(L);
    1862             : }
    1863             : static GEN
    1864         105 : lfununextvec(GEN x)
    1865             : {
    1866         105 :   if (typ(x)==t_SER) return lfununext(x);
    1867           0 :   pari_APPLY_same(lfununext(gel(x,i)));
    1868             : }
    1869             : 
    1870             : /* assume lmisc is an linit, s went through get_domain and s/bitprec belong
    1871             :  * to domain */
    1872             : static GEN
    1873        9856 : lfun_OK(GEN linit, GEN s, GEN sdom, long bitprec)
    1874             : {
    1875        9856 :   GEN N, gas, S, FVga, res, ss = s;
    1876        9856 :   long prec = nbits2prec(bitprec), ext;
    1877             : 
    1878        9856 :   FVga = lfun_get_factgammavec(linit_get_tech(linit));
    1879        9856 :   S = lfunlambda_OK(linit, s, sdom, bitprec);
    1880        9856 :   if (is_ser(S))
    1881             :   {
    1882        1645 :     GEN r = typ(S)==t_SER ? S : gel(S,1);
    1883        1645 :     long d = lg(r) - 2 + fracgammadegree(FVga);
    1884        1645 :     if (typ(s) == t_SER)
    1885        1316 :       ss = sertoser(s, d);
    1886             :     else
    1887         329 :       ss = deg1ser_shallow(gen_1, s, varn(r), d);
    1888             :   }
    1889        9856 :   gas = gammafactproduct(FVga, ss, &ext, prec);
    1890        9856 :   N = ldata_get_conductor(linit_get_ldata(linit));
    1891        9856 :   res = gdiv(S, gmul(gpow(N, gdivgu(ss, 2), prec), gas));
    1892        9856 :   if (typ(s) != t_SER && is_ser(res)) res = lfunservec(res);
    1893        9485 :   else if (ext) res = lfununextvec(res);
    1894        9856 :   return gprec_w(res, prec);
    1895             : }
    1896             : 
    1897             : GEN
    1898       13139 : lfun(GEN lmisc, GEN s, long bitprec)
    1899             : {
    1900       13139 :   pari_sp av = avma;
    1901       13139 :   GEN linit = NULL, ldata, dom, z;
    1902             :   long der;
    1903       13139 :   s = get_domain(s, &dom, &der);
    1904       13139 :   if (der && typ(s) != t_SER)
    1905             :   {
    1906           0 :     if (lfunspec_OK(lmisc, s, &ldata))
    1907             :     {
    1908           0 :       linit = lfuninit(lmisc, cgetg(1,t_VEC), 0, bitprec);
    1909           0 :       return derivnumk((void*)linit, (GEN(*)(void*,GEN,long))&lfun,
    1910             :                        s, stoi(der), nbits2prec(bitprec));
    1911             :     }
    1912             :   }
    1913             :   else
    1914             :   {
    1915       13139 :     long t = lfunspec_OK(lmisc, s, &ldata);
    1916       13132 :     if (t == 1)
    1917             :     { /* special value ? */
    1918        3353 :       long D = itos_or_0(gel(ldata_get_an(ldata), 2)), ss = itos(s);
    1919        3353 :       if (D)
    1920             :       {
    1921        3353 :         if (ss <= 0) return lfunquadneg(D, ss);
    1922             :         /* ss > 0 */
    1923         770 :         if ((!odd(ss) && D > 0) || (odd(ss) && D < 0))
    1924             :         {
    1925         672 :           long prec = nbits2prec(bitprec), q = labs(D);
    1926         672 :           ss = 1 - ss; /* <= 0 */
    1927         672 :           z = powrs(divrs(mppi(prec + EXTRAPREC64), q), 1-ss);
    1928         672 :           z = mulrr(shiftr(z, -ss), sqrtr_abs(utor(q, prec)));
    1929         672 :           z = gdiv(z, mpfactr(-ss, prec));
    1930         672 :           if (smodss(ss, 4) > 1) togglesign(z);
    1931         672 :           return gmul(z, lfunquadneg(D, ss));
    1932             :         }
    1933             :       }
    1934             :     }
    1935        9877 :     if (is_linit(lmisc)) linit = lmisc; else lmisc = ldata;
    1936        9877 :     if (t == 2)
    1937        1211 :       return gc_GEN(av, linit? _product(&lfun, linit, s, bitprec)
    1938         231 :                                    : lfunlarge(ldata, s, bitprec));
    1939             :   }
    1940        8897 :   linit = lfuninit(lmisc, dom, der, bitprec);
    1941        8883 :   z = lfun_OK(linit, s, dom, bitprec);
    1942        8883 :   return gc_GEN(av, z);
    1943             : }
    1944             : 
    1945             : /* given a t_SER a+x*s(x), return x*s(x), shallow */
    1946             : static GEN
    1947          42 : sersplit1(GEN s, GEN *head)
    1948             : {
    1949          42 :   long i, l = lg(s);
    1950             :   GEN y;
    1951          42 :   *head = simplify_shallow(mysercoeff(s, 0));
    1952          42 :   if (valser(s) > 0) return s;
    1953          28 :   y = cgetg(l-1, t_SER); y[1] = s[1];
    1954          28 :   setvalser(y, 1);
    1955         140 :   for (i=3; i < l; i++) gel(y,i-1) = gel(s,i);
    1956          28 :   return normalizeser(y);
    1957             : }
    1958             : 
    1959             : /* order of pole of Lambda at s (0 if regular point) */
    1960             : static long
    1961        2226 : lfunlambdaord(GEN linit, GEN s)
    1962             : {
    1963        2226 :   GEN tech = linit_get_tech(linit);
    1964        2226 :   if (linit_get_type(linit)==t_LDESC_PRODUCT)
    1965             :   {
    1966         287 :     GEN v = lfunprod_get_fact(linit_get_tech(linit));
    1967         287 :     GEN F = gel(v, 1), E = gel(v, 2), C = gel(v, 3);
    1968         287 :     long i, ex = 0, l = lg(F);
    1969         980 :     for (i = 1; i < l; i++)
    1970         693 :       ex += lfunlambdaord(gel(F,i), s) * (E[i]+C[i]);
    1971         287 :     return ex;
    1972             :   }
    1973        1939 :   if (ldata_get_residue(linit_get_ldata(linit)))
    1974             :   {
    1975         595 :     GEN r = lfunpoleresidue(lfun_get_Residue(tech), s);
    1976         595 :     if (r) return lg(r)-2;
    1977             :   }
    1978        1778 :   return 0;
    1979             : }
    1980             : 
    1981             : static GEN
    1982         126 : derser(GEN res, long m)
    1983             : {
    1984         126 :   long v = valser(res);
    1985         126 :   if (v > m) return gen_0;
    1986         126 :   if (v >= 0)
    1987         126 :     return gmul(mysercoeff(res, m), mpfact(m));
    1988             :   else
    1989           0 :     return derivn(res, m, -1);
    1990             : }
    1991             : 
    1992             : static GEN
    1993         189 : derservec(GEN x, long m) { pari_APPLY_same(derser(gel(x,i),m)) }
    1994             : 
    1995             : /* derivative of order m > 0 of L (flag = 0) or Lambda (flag = 1) */
    1996             : static GEN
    1997        1624 : lfunderiv(GEN lmisc, long m, GEN s, long flag, long bitprec)
    1998             : {
    1999        1624 :   pari_sp ltop = avma;
    2000        1624 :   GEN res, S = NULL, linit, ldata, dom;
    2001        1624 :   long der, prec = nbits2prec(bitprec);
    2002        1624 :   if (m <= 0) pari_err_DOMAIN("lfun", "D", "<=", gen_0, stoi(m));
    2003        1617 :   s = get_domain(s, &dom, &der);
    2004        1617 :   if (typ(s) != t_SER && lfunspec_OK(lmisc, s, &ldata) == 2)
    2005             :   {
    2006          28 :     linit = lfuninit(lmisc, cgetg(1,t_VEC), 0, bitprec);
    2007          28 :     return derivnumk((void*)linit, (GEN(*)(void*,GEN,long))&lfun,
    2008             :                      s, stoi(der + m), prec);
    2009             :   }
    2010        1589 :   linit = lfuninit(lmisc, dom, der + m, bitprec);
    2011        1589 :   if (lg(lfun_get_dom(linit_get_tech(linit))) == 1)
    2012          14 :     pari_err_IMPL("domain = [] for derivatives in lfuninit");
    2013        1575 :   if (typ(s) == t_SER)
    2014             :   {
    2015             :     GEN a;
    2016          42 :     if (valser(s) < 0) pari_err_DOMAIN("lfun","valuation", "<", gen_0, s);
    2017          42 :     S = sersplit1(s, &a);
    2018          42 :     s = deg1ser_shallow(gen_1, a, varn(S), m + ceildivuu(lg(s)-2, valser(S)));
    2019             :   }
    2020             :   else
    2021             :   {
    2022        1533 :     long e = lfunlambdaord(linit, s) + m + 1;
    2023             :     /* HACK: pretend lfuninit was done to right accuracy */
    2024        1533 :     if (gequal0(s)) { s = gen_0; e--; }
    2025        1533 :     s = deg1ser_shallow(gen_1, s, 0, e);
    2026             :   }
    2027        1575 :   res = flag ? lfunlambda_OK(linit, s, dom, bitprec):
    2028         973 :                lfun_OK(linit, s, dom, bitprec);
    2029        1575 :   if (S)
    2030          42 :     res = gsubst(derivn(res, m, -1), varn(S), S);
    2031        1533 :   else if (typ(res)==t_SER)
    2032             :   {
    2033        1470 :     long v = valser(res);
    2034        1470 :     if (v > m) { set_avma(ltop); return gen_0; }
    2035        1456 :     if (v >= 0)
    2036        1330 :       res = gmul(mysercoeff(res, m), mpfact(m));
    2037             :     else
    2038         126 :       res = derivn(res, m, -1);
    2039             :   }
    2040          63 :   else if (is_ser(res))
    2041          63 :     res = derservec(res, m);
    2042        1561 :   return gc_GEN(ltop, gprec_w(res, prec));
    2043             : }
    2044             : 
    2045             : GEN
    2046        1512 : lfunlambda0(GEN lmisc, GEN s, long der, long bitprec)
    2047             : {
    2048         623 :   return der? lfunderiv(lmisc, der, s, 1, bitprec)
    2049        2128 :             : lfunlambda(lmisc, s, bitprec);
    2050             : }
    2051             : 
    2052             : GEN
    2053        6874 : lfun0(GEN lmisc, GEN s, long der, long bitprec)
    2054             : {
    2055        1001 :   return der? lfunderiv(lmisc, der, s, 0, bitprec)
    2056        7861 :             : lfun(lmisc, s, bitprec);
    2057             : }
    2058             : 
    2059             : GEN
    2060       19354 : lfunhardy(GEN lmisc, GEN t, long bitprec)
    2061             : {
    2062       19354 :   pari_sp ltop = avma;
    2063       19354 :   long prec = nbits2prec(bitprec), d, isbig = 0;
    2064             :   GEN linit, h, ldata, tech, w2, k2, E, a, argz, z;
    2065             : 
    2066       19354 :   switch(typ(t))
    2067             :   {
    2068       19347 :     case t_INT: case t_FRAC: case t_REAL: break;
    2069           7 :     default: pari_err_TYPE("lfunhardy",t);
    2070             :   }
    2071       19347 :   if (lfunspec_OK(lmisc, mkcomplex(gen_0, t), &ldata) == 2)
    2072             :   {
    2073         868 :     long B = bitprec + maxss(gexpo(t), 0);
    2074         868 :     GEN L = NULL;
    2075         868 :     isbig = 1;
    2076         868 :     k2 = ghalf;
    2077         868 :     z = mkcomplex(k2, t);
    2078         868 :     if (is_linit(lmisc))
    2079             :     {
    2080         742 :       linit = lmisc;
    2081         742 :       if (linit_get_type(linit) == t_LDESC_PRODUCT)
    2082          14 :         L = mkvec(linit);/*HACK*/
    2083             :     }
    2084             :     else
    2085             :     {
    2086         126 :       linit = lfunnoinit(ldata, B);
    2087         126 :       ldata = linit_get_ldata(linit); /* make sure eno is included */
    2088             :     }
    2089         868 :     h = lfunloglambdalarge(L? L: ldata, gprec_w(z, nbits2prec(B)), B);
    2090         868 :     tech = linit_get_tech(linit);
    2091             :   }
    2092             :   else
    2093             :   {
    2094       18479 :     GEN k = ldata_get_k(ldata);
    2095       18479 :     GEN dom = mkvec3(gmul2n(k, -1), gen_0, gabs(t,LOWDEFAULTPREC));
    2096       18479 :     if (!is_linit(lmisc)) lmisc = ldata;
    2097       18479 :     linit = lfuninit(lmisc, dom, 0, bitprec);
    2098       18479 :     tech = linit_get_tech(linit);
    2099       18479 :     k2 = lfun_get_k2(tech);
    2100       18479 :     z = mkcomplex(k2, t);
    2101       18479 :     h = lfunlambda_OK(linit, z, dom, bitprec);
    2102             :   }
    2103       19347 :   w2 = lfun_get_w2(tech);
    2104       19347 :   E = lfun_get_expot(tech); /* 4E = d(k2 - 1) + real(vecsum(Vga)) */
    2105       19347 :   d = ldata_get_degree(ldata);
    2106             :   /* more accurate than garg: k/2 in Q */
    2107       19347 :   argz = gequal0(k2)? Pi2n(-1, prec): gatan(gdiv(t, k2), prec);
    2108       19347 :   prec = precision(argz);
    2109             :   /* prec may have increased: don't lose accuracy if |z|^2 is exact */
    2110       19347 :   a = gsub(gmulsg(d, gmul(t, gmul2n(argz,-1))),
    2111             :            gmul(E, glog(gnorm(z),prec)));
    2112       19347 :   if (!isint1(w2) && typ(ldata_get_dual(ldata))==t_INT)
    2113       16205 :     h = isbig ? gadd(h, glog(w2, prec)) : mulrealvec(h, w2);
    2114       19347 :   if (typ(h) == t_COMPLEX && gexpo(imag_i(h)) < -(bitprec >> 1))
    2115        2533 :     h = real_i(h);
    2116       19347 :   if (isbig) h = greal(gexp(gadd(h, a), prec));
    2117       18479 :   else h = gmul(h, gexp(a, prec));
    2118       19347 :   return gc_upto(ltop, h);
    2119             : }
    2120             : 
    2121             : /* L = log(t); return  \sum_{i = 0}^{v-1}  R[-i-1] L^i/i! */
    2122             : static GEN
    2123        1918 : theta_pole_contrib(GEN R, long v, GEN L)
    2124             : {
    2125        1918 :   GEN s = mysercoeff(R,-v);
    2126             :   long i;
    2127        2023 :   for (i = v-1; i >= 1; i--)
    2128         105 :     s = gadd(mysercoeff(R,-i), gdivgu(gmul(s,L), i));
    2129        1918 :   return s;
    2130             : }
    2131             : /* subtract successively rather than adding everything then subtracting.
    2132             :  * The polar part is "large" and suffers from cancellation: a little stabler
    2133             :  * this way */
    2134             : static GEN
    2135        6951 : theta_add_polar_part(GEN S, GEN R, GEN t, long prec)
    2136             : {
    2137        6951 :   GEN logt = NULL;
    2138        6951 :   long j, l = lg(R);
    2139        8869 :   for (j = 1; j < l; j++)
    2140             :   {
    2141        1918 :     GEN Rj = gel(R,j), b = gel(Rj,1), Rb = gel(Rj,2);
    2142        1918 :     long v = -valser(Rb);
    2143        1918 :     if (v > 1 && !logt) logt = glog(t, prec);
    2144        1918 :     S = gsub(S, gmul(theta_pole_contrib(Rb,v,logt), gpow(t,b,prec)));
    2145             :   }
    2146        6951 :   return S;
    2147             : }
    2148             : 
    2149             : static long
    2150        3598 : lfuncheckfeq_i(GEN theta, GEN thetad, GEN t0, GEN t0i, long bitprec)
    2151             : {
    2152        3598 :   GEN ldata = linit_get_ldata(theta);
    2153             :   GEN S0, S0i, w, eno;
    2154        3598 :   long prec = nbits2prec(bitprec);
    2155        3598 :   if (thetad)
    2156          70 :     S0 = lfuntheta(thetad, t0, 0, bitprec);
    2157             :   else
    2158        3528 :     S0 = conj_i(lfuntheta(theta, conj_i(t0), 0, bitprec));
    2159        3598 :   S0i = lfuntheta(theta, t0i, 0, bitprec);
    2160             : 
    2161        3598 :   eno = ldata_get_rootno(ldata);
    2162        3598 :   if (ldata_get_residue(ldata))
    2163             :   {
    2164         959 :     GEN R = theta_get_R(linit_get_tech(theta));
    2165         959 :     if (gequal0(R))
    2166             :     {
    2167             :       GEN v, r;
    2168         105 :       long t = ldata_get_type(ldata);
    2169         105 :       if (t == t_LFUN_NF || t == t_LFUN_ABELREL)
    2170             :       { /* inefficient since theta not needed; no need to optimize for this
    2171             :            (artificial) query [e.g. lfuncheckfeq(t_POL)] */
    2172          42 :         GEN L = lfuninit(ldata,zerovec(3),0,bitprec);
    2173          42 :         return lfuncheckfeq(L,t0,bitprec);
    2174             :       }
    2175          63 :       v = lfunrootres(theta, bitprec);
    2176          63 :       r = gel(v,1);
    2177          63 :       if (gequal0(eno)) eno = gel(v,3);
    2178          63 :       R = lfunrtoR_i(ldata, r, eno, nbits2prec(bitprec));
    2179             :     }
    2180         917 :     S0i = theta_add_polar_part(S0i, R, t0, prec);
    2181             :   }
    2182        3556 :   if (gequal0(S0i) || gequal0(S0)) pari_err_PREC("lfuncheckfeq");
    2183             : 
    2184        3556 :   w = gdivvec(S0i, gmul(S0, gpow(t0, ldata_get_k(ldata), prec)));
    2185             :   /* missing rootno: guess it */
    2186        3556 :   if (gequal0(eno)) eno = lfunrootno(theta, bitprec);
    2187        3556 :   w = gsubvec(w, eno);
    2188        3556 :   if (thetad) w = gdivvec(w, eno); /* |eno| may be large in non-dual case */
    2189        3556 :   return gexpo(w);
    2190             : }
    2191             : 
    2192             : /* Check whether the coefficients, conductor, weight, polar part and root
    2193             :  * number are compatible with the functional equation at t0 and 1/t0.
    2194             :  * Different from lfunrootres. */
    2195             : long
    2196        3731 : lfuncheckfeq(GEN lmisc, GEN t0, long bitprec)
    2197             : {
    2198             :   GEN ldata, theta, thetad, t0i;
    2199             :   pari_sp av;
    2200             : 
    2201        3731 :   if (is_linit(lmisc) && linit_get_type(lmisc)==t_LDESC_PRODUCT)
    2202             :   {
    2203         168 :     GEN v = lfunprod_get_fact(linit_get_tech(lmisc)), F = gel(v,1);
    2204         168 :     long i, b = -bitprec, l = lg(F);
    2205         560 :     for (i = 1; i < l; i++) b = maxss(b, lfuncheckfeq(gel(F,i), t0, bitprec));
    2206         168 :     return b;
    2207             :   }
    2208        3563 :   av = avma;
    2209        3563 :   if (!t0)
    2210             :   { /* ~Pi/3 + I/7, some random complex number */
    2211        3388 :     t0 = mkcomplex(uutoQ(355,339), uutoQ(1,7));
    2212        3388 :     t0i = ginv(t0);
    2213             :   }
    2214         175 :   else if (gcmpgs(gnorm(t0), 1) < 0) { t0i = t0; t0 = ginv(t0); }
    2215         119 :   else t0i = ginv(t0);
    2216             :   /* |t0| >= 1 */
    2217        3563 :   theta = lfunthetacheckinit(lmisc, t0i, 0, bitprec);
    2218        3556 :   ldata = linit_get_ldata(theta);
    2219        3556 :   thetad = theta_dual(theta, ldata_get_dual(ldata));
    2220        3556 :   return gc_long(av, lfuncheckfeq_i(theta, thetad, t0, t0i, bitprec));
    2221             : }
    2222             : 
    2223             : /*******************************************************************/
    2224             : /*       Compute root number and residues                          */
    2225             : /*******************************************************************/
    2226             : /* round root number to \pm 1 if close to integer. */
    2227             : static GEN
    2228        6384 : ropm1(GEN w, long prec)
    2229             : {
    2230             :   long e;
    2231             :   GEN r;
    2232        6384 :   if (typ(w) == t_INT) return w;
    2233        5978 :   r = grndtoi(w, &e);
    2234        5978 :   return (e < -prec/2)? r: w;
    2235             : }
    2236             : 
    2237             : /* theta for t=1/sqrt(2) and t2==2t simultaneously, saving 25% of the work.
    2238             :  * Assume correct initialization (no thetacheck) */
    2239             : static void
    2240         434 : lfunthetaspec(GEN linit, long bitprec, GEN *pv, GEN *pv2)
    2241             : {
    2242         434 :   pari_sp av = avma, av2;
    2243             :   GEN t, Vga, an, K, ldata, thetainit, v, v2, vroots;
    2244             :   long L, prec, n, d;
    2245             : 
    2246         434 :   ldata = linit_get_ldata(linit);
    2247         434 :   thetainit = linit_get_tech(linit);
    2248         434 :   prec = nbits2prec(bitprec);
    2249         434 :   Vga = ldata_get_gammavec(ldata); d = lg(Vga)-1;
    2250         434 :   if (Vgaeasytheta(Vga))
    2251             :   {
    2252         224 :     GEN v2 = sqrtr(real2n(1, nbits2prec(bitprec)));
    2253         224 :     GEN v = shiftr(v2,-1);
    2254         224 :     *pv = lfuntheta(linit, v,  0, bitprec);
    2255         224 :     *pv2= lfuntheta(linit, v2, 0, bitprec);
    2256         224 :     return;
    2257             :   }
    2258         210 :   an = RgV_kill0( theta_get_an(thetainit) );
    2259         210 :   L = lg(an)-1;
    2260             :   /* to compute theta(1/sqrt(2)) */
    2261         210 :   t = ginv(gsqrt(gmul2n(ldata_get_conductor(ldata), 1), prec));
    2262             :   /* t = 1/sqrt(2N) */
    2263             : 
    2264             :   /* From then on, the code is generic and could be used to compute
    2265             :    * theta(t) / theta(2t) without assuming t = 1/sqrt(2) */
    2266         210 :   K = theta_get_K(thetainit);
    2267         210 :   vroots = mkvroots(d, L, prec);
    2268         210 :   t = gpow(t, gdivgu(gen_2, d), prec); /* rt variant: t->t^(2/d) */
    2269             :   /* v = \sum_{n <= L, n odd} a_n K(nt) */
    2270     1815212 :   for (v = gen_0, n = 1; n <= L; n+=2)
    2271             :   {
    2272     1815002 :     GEN tn, Kn, a = gel(an, n);
    2273             : 
    2274     1815002 :     if (!a) continue;
    2275      113729 :     av2 = avma;
    2276      113729 :     tn = gmul(t, gel(vroots,n));
    2277      113729 :     Kn = gammamellininvrt(K, tn, bitprec);
    2278      113729 :     v = gc_upto(av2, gadd(v, gmul(a,Kn)));
    2279             :   }
    2280             :   /* v += \sum_{n <= L, n even} a_n K(nt), v2 = \sum_{n <= L/2} a_n K(2n t) */
    2281     1815114 :   for (v2 = gen_0, n = 1; n <= L/2; n++)
    2282             :   {
    2283     1814904 :     GEN t2n, K2n, a = gel(an, n), a2 = gel(an,2*n);
    2284             : 
    2285     1814904 :     if (!a && !a2) continue;
    2286      120484 :     av2 = avma;
    2287      120484 :     t2n = gmul(t, gel(vroots,2*n));
    2288      120484 :     K2n = gc_upto(av2, gammamellininvrt(K, t2n, bitprec));
    2289      120484 :     if (a) v2 = gadd(v2, gmul(a, K2n));
    2290      120484 :     if (a2) v = gadd(v,  gmul(a2,K2n));
    2291             :   }
    2292         210 :   *pv = v;
    2293         210 :   *pv2 = v2;
    2294         210 :   (void)gc_all(av, 2, pv,pv2);
    2295             : }
    2296             : 
    2297             : static GEN
    2298         413 : Rtor(GEN a, GEN R, GEN ldata, long prec)
    2299             : {
    2300         413 :   GEN FVga = gammafactor(ldata_get_gammavec(ldata));
    2301         413 :   GEN Na = gpow(ldata_get_conductor(ldata), gdivgu(a,2), prec);
    2302             :   long ext;
    2303         413 :   return gdiv(R, gmul(Na, gammafactproduct(FVga, a, &ext, prec)));
    2304             : }
    2305             : 
    2306             : /* v = theta~(t), vi = theta(1/t) */
    2307             : static GEN
    2308        6034 : get_eno(GEN R, GEN k, GEN t, GEN v, GEN vi, long vx, long bitprec, long force)
    2309             : {
    2310        6034 :   long prec = nbits2prec(bitprec);
    2311        6034 :   GEN a0, a1, S = deg1pol(gmul(gpow(t,k,prec), gneg(v)), vi, vx);
    2312             : 
    2313        6034 :   S = theta_add_polar_part(S, R, t, prec);
    2314        6034 :   if (typ(S) != t_POL || degpol(S) != 1) return NULL;
    2315        6034 :   a1 = gel(S,3); if (!force && gexpo(a1) < -bitprec/4) return NULL;
    2316        5971 :   a0 = gel(S,2);
    2317        5971 :   return gdivvec(a0, gneg(a1));
    2318             : 
    2319             : }
    2320             : /* Return w using theta(1/t) - w t^k \bar{theta}(t) = polar_part(t,w).
    2321             :  * The full Taylor expansion of L must be known */
    2322             : GEN
    2323        5971 : lfunrootno(GEN linit, long bitprec)
    2324             : {
    2325             :   GEN ldata, t, eno, v, vi, R, thetad;
    2326        5971 :   long c = 0, prec = nbits2prec(bitprec), vx = fetch_var();
    2327             :   GEN k;
    2328             :   pari_sp av;
    2329             : 
    2330             :   /* initialize for t > 1/sqrt(2) */
    2331        5971 :   linit = lfunthetacheckinit(linit, dbltor(sqrt(0.5)), 0, bitprec);
    2332        5971 :   ldata = linit_get_ldata(linit);
    2333        5971 :   k = ldata_get_k(ldata);
    2334        5985 :   R = ldata_get_residue(ldata)? lfunrtoR_eno(ldata, pol_x(vx), prec)
    2335        5971 :                               : cgetg(1, t_VEC);
    2336        5971 :   t = gen_1;
    2337        5971 :   v = lfuntheta(linit, t, 0, bitprec);
    2338        5971 :   thetad = theta_dual(linit, ldata_get_dual(ldata));
    2339        5971 :   vi = !thetad ? conj_i(v): lfuntheta(thetad, t, 0, bitprec);
    2340        5971 :   eno = get_eno(R,k,t,vi,v, vx, bitprec, 0);
    2341        5971 :   if (!eno && !thetad)
    2342             :   { /* t = sqrt(2), vi = theta(1/t), v = theta(t) */
    2343          21 :     lfunthetaspec(linit, bitprec, &vi, &v);
    2344          21 :     t = sqrtr(utor(2, prec));
    2345          21 :     eno = get_eno(R,k,t,conj_i(v),vi, vx, bitprec, 0);
    2346             :   }
    2347        5971 :   av = avma;
    2348        6013 :   while (!eno)
    2349             :   {
    2350          42 :     t = addsr(1, shiftr(utor(pari_rand(), prec), -2-BITS_IN_LONG));
    2351             :     /* t in [1,1.25[ */
    2352           0 :     v = thetad? lfuntheta(thetad, t, 0, bitprec)
    2353          42 :               : conj_i(lfuntheta(linit, t, 0, bitprec));
    2354          42 :     vi = lfuntheta(linit, ginv(t), 0, bitprec);
    2355          42 :     eno = get_eno(R,k,t,v,vi, vx, bitprec, c++ == 5);
    2356          42 :     set_avma(av);
    2357             :   }
    2358        5971 :   delete_var(); return ropm1(eno,prec);
    2359             : }
    2360             : 
    2361             : /* Find root number and/or residues when L-function coefficients and
    2362             :    conductor are known. For the moment at most a single residue allowed. */
    2363             : GEN
    2364        6804 : lfunrootres(GEN data, long bitprec)
    2365             : {
    2366        6804 :   pari_sp ltop = avma;
    2367             :   GEN k, w, r, R, a, b, e, v, v2, be, ldata, linit;
    2368             :   long prec;
    2369             : 
    2370        6804 :   ldata = lfunmisc_to_ldata_shallow(data);
    2371        6804 :   r = ldata_get_residue(ldata);
    2372        6804 :   k = ldata_get_k(ldata);
    2373        6804 :   w = ldata_get_rootno(ldata);
    2374        6804 :   if (r) r = normalize_simple_pole(r, k);
    2375        6804 :   if (!r || residues_known(r))
    2376             :   {
    2377        6391 :     if (isintzero(w)) w = lfunrootno(data, bitprec);
    2378        6391 :     if (!r)
    2379        4284 :       r = R = gen_0;
    2380             :     else
    2381        2107 :       R = lfunrtoR_eno(ldata, w, nbits2prec(bitprec));
    2382        6391 :     return gc_GEN(ltop, mkvec3(r, R, w));
    2383             :   }
    2384         413 :   linit = lfunthetacheckinit(data, dbltor(sqrt(0.5)), 0, bitprec);
    2385         413 :   prec = nbits2prec(bitprec);
    2386         413 :   if (lg(r) > 2) pari_err_IMPL("multiple poles in lfunrootres");
    2387             :   /* Now residue unknown, and r = [[be,0]]. */
    2388         413 :   be = gmael(r, 1, 1);
    2389         413 :   if (ldata_isreal(ldata) && gequalm1(w))
    2390           0 :     R = lfuntheta(linit, gen_1, 0, bitprec);
    2391             :   else
    2392             :   {
    2393         413 :     GEN p2k = gpow(gen_2,k,prec);
    2394         413 :     lfunthetaspec(linit, bitprec, &v2, &v);
    2395         413 :     if (gequal(gmulsg(2, be), k)) pari_err_IMPL("pole at k/2 in lfunrootres");
    2396         413 :     if (gequal(be, k))
    2397             :     {
    2398         147 :       a = conj_i(gsub(gmul(p2k, v), v2));
    2399         147 :       b = subiu(p2k, 1);
    2400         147 :       e = gmul(gsqrt(p2k, prec), gsub(v2, v));
    2401             :     }
    2402             :     else
    2403             :     {
    2404         266 :       GEN tk2 = gsqrt(p2k, prec);
    2405         266 :       GEN tbe = gpow(gen_2, be, prec);
    2406         266 :       GEN tkbe = gpow(gen_2, gdivgu(gsub(k, be), 2), prec);
    2407         266 :       a = conj_i(gsub(gmul(tbe, v), v2));
    2408         266 :       b = gsub(gdiv(tbe, tkbe), tkbe);
    2409         266 :       e = gsub(gmul(gdiv(tbe, tk2), v2), gmul(tk2, v));
    2410             :     }
    2411         413 :     if (isintzero(w))
    2412             :     { /* Now residue unknown, r = [[be,0]], and w unknown. */
    2413           7 :       GEN t0  = mkfrac(utoi(11),utoi(10));
    2414           7 :       GEN th1 = lfuntheta(linit, t0,  0, bitprec);
    2415           7 :       GEN th2 = lfuntheta(linit, ginv(t0), 0, bitprec);
    2416           7 :       GEN tbe = gpow(t0, gmulsg(2, be), prec);
    2417           7 :       GEN tkbe = gpow(t0, gsub(k, be), prec);
    2418           7 :       GEN tk2 = gpow(t0, k, prec);
    2419           7 :       GEN c = conj_i(gsub(gmul(tbe, th1), th2));
    2420           7 :       GEN d = gsub(gdiv(tbe, tkbe), tkbe);
    2421           7 :       GEN f = gsub(gmul(gdiv(tbe, tk2), th2), gmul(tk2, th1));
    2422           7 :       GEN D = gsub(gmul(a, d), gmul(b, c));
    2423           7 :       w = gdiv(gsub(gmul(d, e), gmul(b, f)), D);
    2424             :     }
    2425         413 :     w = ropm1(w, prec);
    2426         413 :     R = gdiv(gsub(e, gmul(a, w)), b);
    2427             :   }
    2428         413 :   r = normalize_simple_pole(Rtor(be, R, ldata, prec), be);
    2429         413 :   R = lfunrtoR_i(ldata, r, w, prec);
    2430         413 :   return gc_GEN(ltop, mkvec3(r, R, w));
    2431             : }
    2432             : 
    2433             : /*******************************************************************/
    2434             : /*                           Zeros                                 */
    2435             : /*******************************************************************/
    2436             : struct lhardyz_t {
    2437             :   long bitprec, prec;
    2438             :   GEN linit;
    2439             : };
    2440             : 
    2441             : static GEN
    2442       18570 : lfunhardyzeros(void *E, GEN t)
    2443             : {
    2444       18570 :   struct lhardyz_t *S = (struct lhardyz_t*)E;
    2445       18570 :   GEN z = gprec_wensure(lfunhardy(S->linit, t, S->bitprec), S->prec);
    2446       18570 :   return typ(z) == t_VEC ? RgV_prod(z): z;
    2447             : }
    2448             : 
    2449             : /* initialize for computation on critical line up to height h, zero
    2450             :  * of order <= m */
    2451             : static GEN
    2452         560 : lfuncenterinit(GEN lmisc, double h, long m, long bitprec)
    2453             : {
    2454         560 :   GEN ldata = lfunmisc_to_ldata_shallow(lmisc);
    2455         560 :   if (m < 0)
    2456             :   { /* choose a sensible default */
    2457         560 :     m = 4;
    2458         560 :     if (is_linit(lmisc) && linit_get_type(lmisc) == t_LDESC_INIT)
    2459             :     {
    2460         469 :       GEN domain = lfun_get_domain(linit_get_tech(lmisc));
    2461         469 :       m = domain_get_der(domain);
    2462             :     }
    2463             :   }
    2464         560 :   if (is_dirichlet(ldata)) m = 0;
    2465         560 :   return lfuninit(lmisc, mkvec(dbltor(h)), m, bitprec);
    2466             : }
    2467             : 
    2468             : long
    2469         553 : lfunorderzero(GEN lmisc, long m, long bitprec)
    2470             : {
    2471         553 :   pari_sp ltop = avma;
    2472             :   GEN eno, ldata, linit, k2;
    2473             :   long G, c0, c, st;
    2474             : 
    2475         553 :   if (is_linit(lmisc) && linit_get_type(lmisc) == t_LDESC_PRODUCT)
    2476             :   {
    2477          84 :     GEN M = gmael(linit_get_tech(lmisc), 2,1);
    2478          84 :     long i, l = lg(M);
    2479         280 :     for (c=0, i=1; i < l; i++) c += lfunorderzero(gel(M,i), m, bitprec);
    2480          84 :     return c;
    2481             :   }
    2482         469 :   linit = lfuncenterinit(lmisc, 0, m, bitprec);
    2483         469 :   ldata = linit_get_ldata(linit);
    2484         469 :   eno = ldata_get_rootno(ldata);
    2485         469 :   k2 = gmul2n(ldata_get_k(ldata), -1);
    2486         469 :   G = -bitprec/2;
    2487         469 :   c0 = 0; st = 1;
    2488         469 :   if (typ(eno) == t_VEC)
    2489             :   {
    2490          42 :     long i, l = lg(eno), cnt = l-1, s = 0;
    2491          42 :     GEN v = zero_zv(l-1);
    2492          42 :     if (ldata_isreal(ldata)) st = 2;
    2493          84 :     for (c = c0; cnt; c += st)
    2494             :     {
    2495          42 :       GEN L = lfun0(linit, k2, c, bitprec);
    2496         154 :       for (i = 1; i < l; i++)
    2497             :       {
    2498         112 :         if (v[i]==0 && gexpo(gel(L,i)) > G)
    2499             :         {
    2500         112 :           v[i] = c; cnt--; s += c;
    2501             :         }
    2502             :       }
    2503             :     }
    2504          42 :     return gc_long(ltop,s);
    2505             :   }
    2506             :   else
    2507             :   {
    2508         427 :     if (ldata_isreal(ldata)) { st = 2; if (!gequal1(eno)) c0 = 1; }
    2509         455 :     for (c = c0;; c += st)
    2510         455 :       if (gexpo(lfun0(linit, k2, c, bitprec)) > G) return gc_long(ltop, c);
    2511             :   }
    2512             : }
    2513             : 
    2514             : /* assume T1 * T2 > 0, T1 <= T2 */
    2515             : static void
    2516          98 : lfunzeros_i(struct lhardyz_t *S, GEN *pw, long *ct, GEN T1, GEN T2, long d,
    2517             :             GEN cN, GEN pi2, GEN pi2div, long precinit, long prec)
    2518             : {
    2519          98 :   GEN T = T1, w = *pw;
    2520          98 :   long W = lg(w)-1, s = gsigne(lfunhardyzeros(S, T1));
    2521             :   for(;;)
    2522         427 :   {
    2523         525 :     pari_sp av = avma;
    2524             :     GEN D, T0, z;
    2525         525 :     D = gcmp(T, pi2) < 0? cN
    2526         525 :                         : gadd(cN, gmulsg(d, glog(gdiv(T, pi2), prec)));
    2527         525 :     D = gdiv(pi2div, gmulsg(d, D));
    2528             :     for(;;)
    2529       13482 :     {
    2530             :       long s0;
    2531       14007 :       T0 = T; T = gadd(T, D);
    2532       14007 :       if (gcmp(T, T2) >= 0) T = T2;
    2533       14007 :       s0 = gsigne(lfunhardyzeros(S, T));
    2534       14007 :       if (s0 != s) { s = s0; break; }
    2535       13580 :       if (T == T2) { setlg(w, *ct); *pw = w; return; }
    2536             :     }
    2537         427 :     z = zbrent(S, lfunhardyzeros, T0, T, prec); /* T <= T2 */
    2538         427 :     (void)gc_all(av, 2, &T, &z);
    2539         427 :     if (*ct > W) { W *= 2; w = vec_lengthen(w, W); }
    2540         427 :     if (typ(z) == t_REAL) z  = rtor(z, precinit);
    2541         427 :     gel(w, (*ct)++) = z;
    2542             :   }
    2543             :   setlg(w, *ct); *pw = w;
    2544             : }
    2545             : GEN
    2546          98 : lfunzeros(GEN ldata, GEN lim, long divz, long bitprec)
    2547             : {
    2548          98 :   pari_sp ltop = avma;
    2549             :   GEN linit, pi2, pi2div, cN, w, T, h1, h2;
    2550          98 :   long i, d, NEWD, c, ct, s1, s2, prec, prec0 = nbits2prec(bitprec);
    2551             :   double maxt;
    2552             :   struct lhardyz_t S;
    2553             : 
    2554          98 :   if (is_linit(ldata) && linit_get_type(ldata) == t_LDESC_PRODUCT)
    2555             :   {
    2556           0 :     GEN M = gmael(linit_get_tech(ldata), 2,1);
    2557           0 :     long l = lg(M);
    2558           0 :     w = cgetg(l, t_VEC);
    2559           0 :     for (i = 1; i < l; i++) gel(w,i) = lfunzeros(gel(M,i), lim, divz, bitprec);
    2560           0 :     return gc_upto(ltop, vecsort0(shallowconcat1(w), NULL, 0));
    2561             :   }
    2562          98 :   if (typ(lim) == t_VEC)
    2563             :   {
    2564             :     double H1, H2;
    2565          63 :     if (lg(lim) != 3 || gcmp(gel(lim, 1), gel(lim, 2)) >= 0)
    2566           7 :       pari_err_TYPE("lfunzeros",lim);
    2567          56 :     h1 = gel(lim, 1); H1 = gtodouble(h1);
    2568          56 :     h2 = gel(lim, 2); H2 = gtodouble(h2);
    2569          56 :     maxt = maxdd(fabs(H1), fabs(H2));
    2570          56 :     if (H1 * H2 > 0)
    2571             :     {
    2572          35 :       GEN LDATA = lfunmisc_to_ldata_shallow(ldata);
    2573          35 :       double m = mindd(fabs(H1), fabs(H2));
    2574          35 :       if (is_dirichlet(LDATA) && m > lfuninit_cutoff(LDATA)) maxt = 0;
    2575             :     }
    2576             :   }
    2577             :   else
    2578             :   {
    2579          35 :     if (gcmp(lim, gen_0) <= 0) pari_err_TYPE("lfunzeros",lim);
    2580          35 :     h1 = gen_0;
    2581          35 :     h2 = lim;
    2582          35 :     maxt = gtodouble(h2);
    2583             :   }
    2584          91 :   S.linit = linit = lfuncenterinit(ldata, maxt, -1, bitprec);
    2585          91 :   S.bitprec = bitprec;
    2586          91 :   S.prec = prec0;
    2587          91 :   ldata = linit_get_ldata(linit);
    2588          91 :   d = ldata_get_degree(ldata);
    2589             : 
    2590          91 :   NEWD = minss((long) ceil(bitprec + (M_PI/(4*M_LN2)) * d * maxt),
    2591             :                lfun_get_bitprec(linit_get_tech(linit)));
    2592          91 :   prec = nbits2prec(NEWD);
    2593          91 :   cN = gdiv(ldata_get_conductor(ldata), gpowgs(Pi2n(-1, prec), d));
    2594          91 :   cN = gexpo(cN) >= 0? gaddsg(d, gmulsg(2, glog(cN, prec))): utoi(d);
    2595          91 :   pi2 = Pi2n(1, prec);
    2596          91 :   pi2div = gdivgu(pi2, labs(divz));
    2597          91 :   s1 = gsigne(h1);
    2598          91 :   s2 = gsigne(h2);
    2599          91 :   w = cgetg(100+1, t_VEC); c = 1; ct = 0; T = NULL;
    2600          91 :   if (s1 <= 0 && s2 >= 0)
    2601             :   {
    2602          56 :     GEN r = ldata_get_residue(ldata);
    2603          56 :     if (!r || gequal0(r))
    2604             :     {
    2605          35 :       ct = lfunorderzero(linit, -1, bitprec);
    2606          35 :       if (ct) T = real2n(-prec / (2*ct), prec);
    2607             :     }
    2608             :   }
    2609          91 :   if (s1 <= 0)
    2610             :   {
    2611          63 :     if (s1 < 0)
    2612          21 :       lfunzeros_i(&S, &w, &c, h1, T? negr(T): h2,
    2613             :                   d, cN, pi2, pi2div, prec0, prec);
    2614          63 :     if (ct)
    2615             :     {
    2616          21 :       long n = lg(w)-1;
    2617          21 :       if (c + ct >= n) w = vec_lengthen(w, n + ct);
    2618          84 :       for (i = 1; i <= ct; i++) gel(w,c++) = gen_0;
    2619             :     }
    2620             :   }
    2621          91 :   if (s2 > 0 && (T || s1 >= 0))
    2622          77 :     lfunzeros_i(&S, &w, &c, T? T: h1, h2, d, cN, pi2, pi2div, prec0, prec);
    2623          91 :   return gc_GEN(ltop, w);
    2624             : }
    2625             : 
    2626             : /*******************************************************************/
    2627             : /*       Guess conductor                                           */
    2628             : /*******************************************************************/
    2629             : struct huntcond_t {
    2630             :   GEN k;
    2631             :   GEN theta, thetad;
    2632             :   GEN *pM, *psqrtM, *pMd, *psqrtMd;
    2633             : };
    2634             : 
    2635             : static void
    2636       11962 : condset(struct huntcond_t *S, GEN M, long prec)
    2637             : {
    2638       11962 :   *(S->pM) = M;
    2639       11962 :   *(S->psqrtM) = gsqrt(ginv(M), prec);
    2640       11962 :   if (S->thetad != S->theta)
    2641             :   {
    2642           0 :     *(S->pMd) = *(S->pM);
    2643           0 :     *(S->psqrtMd) = *(S->psqrtM);
    2644             :   }
    2645       11962 : }
    2646             : 
    2647             : /* M should eventually converge to N, the conductor. L has no pole. */
    2648             : static GEN
    2649        6888 : wrap1(void *E, GEN M)
    2650             : {
    2651        6888 :   struct huntcond_t *S = (struct huntcond_t*)E;
    2652             :   GEN thetainit, tk, p1, p1inv;
    2653        6888 :   GEN t = mkfrac(stoi(11), stoi(10));
    2654             :   long prec, bitprec;
    2655             : 
    2656        6888 :   thetainit = linit_get_tech(S->theta);
    2657        6888 :   bitprec = theta_get_bitprec(thetainit);
    2658        6888 :   prec = nbits2prec(bitprec);
    2659        6888 :   condset(S, M, prec);
    2660        6888 :   tk = gpow(t, S->k, prec);
    2661        6888 :   p1 = lfuntheta(S->thetad, t, 0, bitprec);
    2662        6888 :   p1inv = lfuntheta(S->theta, ginv(t), 0, bitprec);
    2663        6888 :   return glog(gabs(gmul(tk, gdiv(p1, p1inv)), prec), prec);
    2664             : }
    2665             : 
    2666             : /* M should eventually converge to N, the conductor. L has a pole. */
    2667             : static GEN
    2668        5032 : wrap2(void *E, GEN M)
    2669             : {
    2670        5032 :   struct huntcond_t *S = (struct huntcond_t*)E;
    2671             :   GEN t1k, t2k, p1, p1inv, p2, p2inv, thetainit, R;
    2672        5032 :   GEN t1 = mkfrac(stoi(11), stoi(10)), t2 = mkfrac(stoi(13), stoi(11));
    2673             :   GEN t1be, t2be, t1bemk, t2bemk, t1kmbe, t2kmbe;
    2674             :   GEN F11, F12, F21, F22, P1, P2, res;
    2675             :   long prec, bitprec;
    2676        5032 :   GEN k = S->k;
    2677             : 
    2678        5032 :   thetainit = linit_get_tech(S->theta);
    2679        5032 :   bitprec = theta_get_bitprec(thetainit);
    2680        5032 :   prec = nbits2prec(bitprec);
    2681        5032 :   condset(S, M, prec);
    2682             : 
    2683        5032 :   p1 = lfuntheta(S->thetad, t1, 0, bitprec);
    2684        5032 :   p2 = lfuntheta(S->thetad, t2, 0, bitprec);
    2685        5032 :   p1inv = lfuntheta(S->theta, ginv(t1), 0, bitprec);
    2686        5032 :   p2inv = lfuntheta(S->theta, ginv(t2), 0, bitprec);
    2687        5032 :   t1k = gpow(t1, k, prec);
    2688        5032 :   t2k = gpow(t2, k, prec);
    2689        5032 :   R = theta_get_R(thetainit);
    2690        5032 :   if (typ(R) == t_VEC)
    2691             :   {
    2692           0 :     GEN be = gmael(R, 1, 1);
    2693           0 :     t1be = gpow(t1, be, prec); t1bemk = gdiv(gsqr(t1be), t1k);
    2694           0 :     t2be = gpow(t2, be, prec); t2bemk = gdiv(gsqr(t2be), t2k);
    2695           0 :     t1kmbe = gdiv(t1k, t1be);
    2696           0 :     t2kmbe = gdiv(t2k, t2be);
    2697             :   }
    2698             :   else
    2699             :   { /* be = k */
    2700        5032 :     t1be = t1k; t1bemk = t1k; t1kmbe = gen_1;
    2701        5032 :     t2be = t2k; t2bemk = t2k; t2kmbe = gen_1;
    2702             :   }
    2703        5032 :   F11 = conj_i(gsub(gmul(gsqr(t1be), p1), p1inv));
    2704        5032 :   F12 = conj_i(gsub(gmul(gsqr(t2be), p2), p2inv));
    2705        5032 :   F21 = gsub(gmul(t1k, p1), gmul(t1bemk, p1inv));
    2706        5032 :   F22 = gsub(gmul(t2k, p2), gmul(t2bemk, p2inv));
    2707        5032 :   P1 = gsub(gmul(t1bemk, t1be), t1kmbe);
    2708        5032 :   P2 = gsub(gmul(t2bemk, t2be), t2kmbe);
    2709        5032 :   res = gdiv(gsub(gmul(P2,F21), gmul(P1,F22)),
    2710             :              gsub(gmul(P2,F11), gmul(P1,F12)));
    2711        5032 :   return glog(gabs(res, prec), prec);
    2712             : }
    2713             : 
    2714             : /* If flag = 0 (default) return all conductors found as integers. If
    2715             : flag = 1, return the approximations, not the integers. If flag = 2,
    2716             : return all, even nonintegers. */
    2717             : 
    2718             : static GEN
    2719          84 : checkconductor(GEN v, long bit, long flag)
    2720             : {
    2721             :   GEN w;
    2722          84 :   long e, j, k, l = lg(v);
    2723          84 :   if (flag == 2) return v;
    2724          84 :   w = cgetg(l, t_VEC);
    2725         322 :   for (j = k = 1; j < l; j++)
    2726             :   {
    2727         238 :     GEN N = grndtoi(gel(v,j), &e);
    2728         238 :     if (e < -bit) gel(w,k++) = flag ? gel(v,j): N;
    2729             :   }
    2730          84 :   if (k == 2) return gel(w,1);
    2731           7 :   setlg(w,k); return w;
    2732             : }
    2733             : 
    2734             : static GEN
    2735          98 : parse_maxcond(GEN maxN)
    2736             : {
    2737             :   GEN M;
    2738          98 :   if (!maxN)
    2739          49 :     M = utoipos(10000);
    2740          49 :   else if (typ(maxN) == t_VEC)
    2741             :   {
    2742          14 :     if (!RgV_is_ZV(maxN)) pari_err_TYPE("lfunconductor",maxN);
    2743          14 :     return ZV_sort_shallow(maxN);
    2744             :   }
    2745             :   else
    2746          35 :     M = maxN;
    2747          84 :   return (typ(M) == t_INT)? addiu(M, 1): gceil(M);
    2748             : }
    2749             : 
    2750             : GEN
    2751          98 : lfunconductor(GEN data, GEN maxcond, long flag, long bitprec)
    2752             : {
    2753             :   struct huntcond_t S;
    2754          98 :   pari_sp av = avma;
    2755          98 :   GEN ldata = lfunmisc_to_ldata_shallow(data);
    2756          98 :   GEN ld, r, v, theta, thetad, M, tdom, t0 = NULL, t0i = NULL;
    2757             :   GEN (*eval)(void *, GEN);
    2758             :   long prec;
    2759          98 :   M = parse_maxcond(maxcond);
    2760          98 :   r = ldata_get_residue(ldata);
    2761          98 :   if (typ(M) == t_VEC) /* select in list */
    2762             :   {
    2763          14 :     if (lg(M) == 1) retgc_const(av, cgetg(1, t_VEC));
    2764           7 :     eval = NULL; tdom = dbltor(0.7);
    2765             :   }
    2766          84 :   else if (!r) { eval = wrap1; tdom = uutoQ(10,11); }
    2767             :   else
    2768             :   {
    2769          21 :     if (typ(r) == t_VEC && lg(r) > 2)
    2770           0 :       pari_err_IMPL("multiple poles in lfunconductor");
    2771          21 :     eval = wrap2; tdom = uutoQ(11,13);
    2772             :   }
    2773          91 :   if (eval) bitprec += bitprec/2;
    2774          91 :   prec = nbits2prec(bitprec);
    2775          91 :   ld = shallowcopy(ldata);
    2776          91 :   gel(ld, 5) = eval? M: veclast(M);
    2777          91 :   theta = lfunthetainit_i(ld, tdom, 0, bitprec);
    2778          91 :   thetad = theta_dual(theta, ldata_get_dual(ldata));
    2779          91 :   gel(theta,3) = shallowcopy(linit_get_tech(theta));
    2780          91 :   S.k = ldata_get_k(ldata);
    2781          91 :   S.theta = theta;
    2782          91 :   S.thetad = thetad? thetad: theta;
    2783          91 :   S.pM = &gel(linit_get_ldata(theta),5);
    2784          91 :   S.psqrtM = &gel(linit_get_tech(theta),7);
    2785          91 :   if (thetad)
    2786             :   {
    2787           0 :     S.pMd = &gel(linit_get_ldata(thetad),5);
    2788           0 :     S.psqrtMd = &gel(linit_get_tech(thetad),7);
    2789             :   }
    2790          91 :   if (!eval)
    2791             :   {
    2792           7 :     long i, besti = 0, beste = -10, l = lg(M);
    2793           7 :     t0 = uutoQ(11,10); t0i = uutoQ(10,11);
    2794          49 :     for (i = 1; i < l; i++)
    2795             :     {
    2796          42 :       pari_sp av2 = avma;
    2797             :       long e;
    2798          42 :       condset(&S, gel(M,i), prec);
    2799          42 :       e = lfuncheckfeq_i(theta, thetad, t0, t0i, bitprec);
    2800          42 :       set_avma(av2);
    2801          42 :       if (e < beste) { beste = e; besti = i; }
    2802          35 :       else if (e == beste) beste = besti = 0; /* tie: forget */
    2803             :     }
    2804           7 :     if (!besti) retgc_const(av, cgetg(1, t_VEC));
    2805           7 :     return gc_GEN(av, mkvec2(gel(M,besti), stoi(beste)));
    2806             :   }
    2807          84 :   v = solvestep((void*)&S, eval, ghalf, M, gen_2, 14, prec);
    2808          84 :   return gc_GEN(av, checkconductor(v, bitprec/2, flag));
    2809             : }
    2810             : 
    2811             : /* assume chi primitive */
    2812             : static GEN
    2813        2863 : znchargauss_i(GEN G, GEN chi, long bitprec)
    2814             : {
    2815        2863 :   GEN z, q, F = znstar_get_N(G);
    2816             :   long prec;
    2817             : 
    2818        2863 :   if (equali1(F)) return gen_1;
    2819        2863 :   prec = nbits2prec(bitprec);
    2820        2863 :   q = sqrtr_abs(itor(F, prec));
    2821        2863 :   z = lfuntheta(mkvec2(G,chi), gen_1, 0, bitprec);
    2822        2863 :   if (gexpo(z) < 10 - bitprec)
    2823             :   {
    2824          28 :     if (equaliu(F,300))
    2825             :     {
    2826          14 :       GEN z = rootsof1u_cx(25, prec);
    2827          14 :       GEN n = znconreyexp(G, chi);
    2828          14 :       if (equaliu(n, 131)) return gmul(q, gpowgs(z,14));
    2829           7 :       if (equaliu(n, 71)) return gmul(q, gpowgs(z,11));
    2830             :     }
    2831          14 :     if (equaliu(F,600))
    2832             :     {
    2833          14 :       GEN z = rootsof1u_cx(25, prec);
    2834          14 :       GEN n = znconreyexp(G, chi);
    2835          14 :       if (equaliu(n, 491)) return gmul(q, gpowgs(z,7));
    2836           7 :       if (equaliu(n, 11)) return gmul(q, gpowgs(z,18));
    2837             :     }
    2838           0 :     pari_err_BUG("znchargauss [ Theta(chi,1) = 0 ]");
    2839             :   }
    2840        2835 :   z = gmul(gdiv(z, conj_i(z)), q);
    2841        2835 :   if (zncharisodd(G,chi)) z = mulcxI(z);
    2842        2835 :   return z;
    2843             : }
    2844             : static GEN
    2845        2863 : Z_radical(GEN N, long *om)
    2846             : {
    2847        2863 :   GEN P = gel(Z_factor(N), 1);
    2848        2863 :   *om = lg(P)-1; return ZV_prod(P);
    2849             : }
    2850             : GEN
    2851        5516 : znchargauss(GEN G, GEN chi, GEN a, long bitprec)
    2852             : {
    2853             :   GEN v, T, N, F, b0, b1, b2, bF, a1, aF, A, r, GF, tau, B, faB, u, S;
    2854        5516 :   long omb0, prec = nbits2prec(bitprec);
    2855        5516 :   pari_sp av = avma;
    2856             : 
    2857        5516 :   if (typ(chi) != t_COL) chi = znconreylog(G,chi);
    2858        5516 :   T = znchartoprimitive(G, chi);
    2859        5516 :   GF  = gel(T,1);
    2860        5516 :   chi = gel(T,2); /* now primitive */
    2861        5516 :   N = znstar_get_N(G);
    2862        5516 :   F = znstar_get_N(GF);
    2863        5516 :   if (equalii(N,F)) b1 = bF = gen_1;
    2864             :   else
    2865             :   {
    2866         245 :     v = Z_ppio(diviiexact(N,F), F);
    2867         245 :     bF = gel(v,2); /* (N/F, F^oo) */
    2868         245 :     b1 = gel(v,3); /* cofactor */
    2869             :   }
    2870        5516 :   if (!a) a = a1 = aF = gen_1;
    2871             :   else
    2872             :   {
    2873        5467 :     if (typ(a) != t_INT) pari_err_TYPE("znchargauss",a);
    2874        5467 :     a = modii(a, N);
    2875        5467 :     if (!signe(a)) { set_avma(av); return is_pm1(F)? eulerphi(N): gen_0; }
    2876        3031 :     v = Z_ppio(a, F);
    2877        3031 :     aF = gel(v,2);
    2878        3031 :     a1 = gel(v,3);
    2879             :   }
    2880        3080 :   if (!equalii(aF, bF)) { set_avma(av); return gen_0; }
    2881        2863 :   b0 = Z_radical(b1, &omb0);
    2882        2863 :   b2 = diviiexact(b1, b0);
    2883        2863 :   A = dvmdii(a1, b2, &r);
    2884        2863 :   if (r != gen_0) { set_avma(av); return gen_0; }
    2885        2863 :   B = gcdii(A,b0); faB = Z_factor(B); /* squarefree */
    2886        2863 :   S = eulerphi(mkvec2(B,faB));
    2887        2863 :   if (odd(omb0 + lg(gel(faB,1))-1)) S = negi(S); /* moebius(b0/B) * phi(B) */
    2888        2863 :   S = mulii(S, mulii(aF,b2));
    2889        2863 :   tau = znchargauss_i(GF, chi, bitprec);
    2890        2863 :   u = Fp_div(b0, A, F);
    2891        2863 :   if (!equali1(u))
    2892             :   {
    2893           7 :     GEN ord = zncharorder(GF, chi), z = rootsof1_cx(ord, prec);
    2894           7 :     tau = gmul(tau, znchareval(GF, chi, u, mkvec2(z,ord)));
    2895             :   }
    2896        2863 :   return gc_upto(av, gmul(tau, S));
    2897             : }

Generated by: LCOV version 1.16