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 - lfunutils.c (source / functions) Hit Total Coverage
Test: PARI/GP v2.13.0 lcov report (development 25825-edc109b529) Lines: 1596 1721 92.7 %
Date: 2020-09-21 06:08:33 Functions: 149 158 94.3 %
Legend: Lines: hit not hit

          Line data    Source code
       1             : /* Copyright (C) 2015  The PARI group.
       2             : 
       3             : This file is part of the PARI/GP package.
       4             : 
       5             : PARI/GP is free software; you can redistribute it and/or modify it under the
       6             : terms of the GNU General Public License as published by the Free Software
       7             : Foundation. It is distributed in the hope that it will be useful, but WITHOUT
       8             : ANY WARRANTY WHATSOEVER.
       9             : 
      10             : Check the License for details. You should have received a copy of it, along
      11             : with the package; see the file 'COPYING'. If not, write to the Free Software
      12             : Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. */
      13             : 
      14             : /********************************************************************/
      15             : /**                                                                **/
      16             : /**                 L-functions: Applications                      **/
      17             : /**                                                                **/
      18             : /********************************************************************/
      19             : 
      20             : #include "pari.h"
      21             : #include "paripriv.h"
      22             : 
      23             : static GEN
      24       13601 : tag(GEN x, long t) { return mkvec2(mkvecsmall(t), x); }
      25             : 
      26             : /* v a t_VEC of length > 1 */
      27             : static int
      28       56595 : is_tagged(GEN v)
      29             : {
      30       56595 :   GEN T = gel(v,1);
      31       56595 :   return (typ(T)==t_VEC && lg(T)==3 && typ(gel(T,1))==t_VECSMALL);
      32             : }
      33             : static void
      34       56595 : checkldata(GEN ldata)
      35             : {
      36             :   GEN vga, w, N;
      37             : #if 0 /* assumed already checked and true */
      38             :   long l = lg(ldata);
      39             :   if (typ(ldata)!=t_VEC || l < 7 || l > 8 || !is_tagged(ldata))
      40             :     pari_err_TYPE("checkldata", ldata);
      41             : #endif
      42       56595 :   vga = ldata_get_gammavec(ldata);
      43       56595 :   if (typ(vga) != t_VEC) pari_err_TYPE("checkldata [gammavec]",vga);
      44       56595 :   w = gel(ldata, 4); /* FIXME */
      45       56595 :   switch(typ(w))
      46             :   {
      47       54810 :     case t_INT: case t_FRAC: break;
      48        1785 :     case t_VEC: if (lg(w) == 3 && is_rational_t(typ(gel(w,1)))) break;
      49           0 :     default: pari_err_TYPE("checkldata [weight]",w);
      50             :   }
      51       56595 :   N = ldata_get_conductor(ldata);
      52       56595 :   if (typ(N) != t_INT) pari_err_TYPE("checkldata [conductor]",N);
      53       56595 : }
      54             : 
      55             : /* data may be either an object (polynomial, elliptic curve, etc...)
      56             :  * or a description vector [an,sd,Vga,k,conductor,rootno,{poles}]. */
      57             : GEN
      58        1645 : lfuncreate(GEN data)
      59             : {
      60        1645 :   long lx = lg(data);
      61        1645 :   if (typ(data)==t_VEC && (lx == 7 || lx == 8))
      62             :   {
      63         574 :     GEN ldata = gcopy(data);
      64         574 :     if (!is_tagged(data))
      65             :     { /* tag first component as t_LFUN_GENERIC */
      66         434 :       gel(ldata, 1) = tag(gel(ldata,1), t_LFUN_GENERIC);
      67         434 :       if (typ(gel(ldata, 2))!=t_INT)
      68          21 :         gel(ldata, 2) = tag(gel(ldata,2), t_LFUN_GENERIC);
      69             :     }
      70         574 :     checkldata(ldata); return ldata;
      71        1071 :   } else if (typ(data)==t_CLOSURE && closure_arity(data)==0)
      72             :   {
      73          14 :     pari_sp av = avma;
      74          14 :     GEN ldata = lfuncreate(closure_callgen0prec(data, DEFAULTPREC));
      75          14 :     gel(ldata,1) = tag(data, t_LFUN_CLOSURE0);
      76          14 :     return gerepilecopy(av, ldata);
      77             :   }
      78        1057 :   return lfunmisc_to_ldata(data);
      79             : }
      80             : 
      81             : /********************************************************************/
      82             : /**                     Simple constructors                        **/
      83             : /********************************************************************/
      84             : 
      85             : static GEN
      86          98 : vecan_conj(GEN an, long n, long prec)
      87             : {
      88          98 :   GEN p1 = ldata_vecan(gel(an,1), n, prec);
      89          98 :   return typ(p1) == t_VEC? conj_i(p1): p1;
      90             : }
      91             : 
      92             : static GEN
      93         315 : vecan_mul(GEN an, long n, long prec)
      94             : {
      95         315 :   GEN p1 = ldata_vecan(gel(an,1), n, prec);
      96         315 :   GEN p2 = ldata_vecan(gel(an,2), n, prec);
      97         315 :   if (typ(p1) == t_VECSMALL) p1 = vecsmall_to_vec(p1);
      98         315 :   if (typ(p2) == t_VECSMALL) p2 = vecsmall_to_vec(p2);
      99         315 :   return dirmul(p1, p2);
     100             : }
     101             : 
     102             : static GEN
     103          77 : lfunconvol(GEN a1, GEN a2)
     104          77 : { return tag(mkvec2(a1, a2), t_LFUN_MUL); }
     105             : 
     106             : static GEN
     107         637 : vecan_div(GEN an, long n, long prec)
     108             : {
     109         637 :   GEN p1 = ldata_vecan(gel(an,1), n, prec);
     110         637 :   GEN p2 = ldata_vecan(gel(an,2), n, prec);
     111         637 :   if (typ(p1) == t_VECSMALL) p1 = vecsmall_to_vec(p1);
     112         637 :   if (typ(p2) == t_VECSMALL) p2 = vecsmall_to_vec(p2);
     113         637 :   return dirdiv(p1, p2);
     114             : }
     115             : 
     116             : static GEN
     117          56 : lfunconvolinv(GEN a1, GEN a2)
     118          56 : { return tag(mkvec2(a1,a2), t_LFUN_DIV); }
     119             : 
     120             : static GEN
     121          70 : lfunconj(GEN a1)
     122          70 : { return tag(mkvec(a1), t_LFUN_CONJ); }
     123             : 
     124             : static GEN
     125         133 : lfuncombdual(GEN (*fun)(GEN, GEN), GEN ldata1, GEN ldata2)
     126             : {
     127         133 :   GEN a1 = ldata_get_an(ldata1), a2 = ldata_get_an(ldata2);
     128         133 :   GEN b1 = ldata_get_dual(ldata1), b2 = ldata_get_dual(ldata2);
     129         133 :   if (typ(b1)==t_INT && typ(b2)==t_INT)
     130         133 :     return utoi(signe(b1) || signe(b2));
     131             :   else
     132             :   {
     133           0 :     if (typ(b1)==t_INT) b1 = signe(b1) ? lfunconj(a1): a1;
     134           0 :     if (typ(b2)==t_INT) b2 = signe(b2) ? lfunconj(a2): a2;
     135           0 :     return fun(b1, b2);
     136             :   }
     137             : }
     138             : 
     139             : static GEN
     140         581 : vecan_twist(GEN an, long n, long prec)
     141             : {
     142         581 :   GEN p1 = ldata_vecan(gel(an,1), n, prec);
     143         581 :   GEN p2 = ldata_vecan(gel(an,2), n, prec);
     144             :   long i;
     145             :   GEN V;
     146         581 :   if (typ(p1) == t_VECSMALL) p1 = vecsmall_to_vec(p1);
     147         581 :   if (typ(p2) == t_VECSMALL) p2 = vecsmall_to_vec(p2);
     148         581 :   V = cgetg(n+1, t_VEC);
     149      390418 :   for(i = 1; i <= n ; i++)
     150      389837 :     gel(V, i) = gmul(gel(p1, i), gel(p2, i));
     151         581 :   return V;
     152             : }
     153             : 
     154             : static GEN
     155         602 : vecan_shift(GEN an, long n, long prec)
     156             : {
     157         602 :   GEN p1 = ldata_vecan(gel(an,1), n, prec);
     158         602 :   GEN s = gel(an,2);
     159             :   long i;
     160             :   GEN V;
     161         602 :   if (typ(p1) == t_VECSMALL) p1 = vecsmall_to_vec(p1);
     162         602 :   V = cgetg(n+1, t_VEC);
     163         602 :   if (typ(s)==t_INT)
     164             :   {
     165         413 :     if (equali1(s))
     166       33285 :       for(i = 1; i <= n ; i++)
     167             :       {
     168       32879 :         GEN gi = gel(p1, i);
     169       32879 :         gel(V, i) = gequal0(gi)? gi: gmulgs(gi, i);
     170             :       }
     171             :     else
     172          49 :       for(i = 1; i <= n ; i++)
     173             :       {
     174          42 :         GEN gi = gel(p1, i);
     175          42 :         gel(V, i) = gequal0(gi)? gi: gmul(gi, powgi(utoi(i), s));
     176             :       }
     177             :   }
     178             :   else
     179             :   {
     180         189 :     GEN D = dirpowers(n, s, prec);
     181        7028 :     for(i = 1; i <= n ; i++)
     182        6839 :       gel(V, i) = gmul(gel(p1,i), gel(D,i));
     183             :   }
     184         602 :   return V;
     185             : }
     186             : 
     187             : static GEN
     188         385 : deg1ser_shallow(GEN a1, GEN a0, long e)
     189         385 : { return RgX_to_ser(deg1pol_shallow(a1, a0, 0), e+2); }
     190             : static GEN
     191         126 : residuetopoles(GEN r, GEN k)
     192             : {
     193         126 :   if (!is_vec_t(typ(r))) r = mkvec(mkvec2(k, r));
     194         126 :   return lfunrtopoles(r);
     195             : }
     196             : static GEN
     197          77 : lfunmulpoles(GEN ldata1, GEN ldata2, long bitprec)
     198             : {
     199          77 :   GEN k = ldata_get_k(ldata1);
     200          77 :   GEN r1 = ldata_get_residue(ldata1);
     201          77 :   GEN r2 = ldata_get_residue(ldata2), r;
     202             :   long l, j;
     203             : 
     204          77 :   if (!r1 && !r2) return NULL;
     205          70 :   if (r1) r1 = residuetopoles(r1, k);
     206          70 :   if (r2) r2 = residuetopoles(r2, k);
     207          70 :   r = r1? (r2? setunion_i(r1, r2): r1): r2;
     208          70 :   l = lg(r); if (l == 1) return NULL;
     209         168 :   for (j = 1; j < l; j++)
     210             :   {
     211          98 :     GEN be = gel(r,j), bx = deg1ser_shallow(gen_1, be, 2);
     212          98 :     GEN z1 = lfun(ldata1,bx,bitprec), z2 = lfun(ldata2,bx,bitprec);
     213          98 :     long e = valp(z1) + valp(z2);
     214          98 :     GEN b = deg1ser_shallow(gen_1, be, 2-e);
     215          98 :     z1 = lfun(ldata1,b,bitprec);
     216          98 :     z2 = lfun(ldata2,b,bitprec);
     217          98 :     gel(r,j) = mkvec2(be, gmul(z1, z2));
     218             :   }
     219          70 :   return r;
     220             : }
     221             : 
     222             : static GEN
     223          77 : lfunmul_k(GEN ldata1, GEN ldata2, GEN k, long bitprec)
     224             : {
     225             :   GEN r, N, Vga, eno, a1a2, b1b2;
     226          77 :   r = lfunmulpoles(ldata1, ldata2, bitprec);
     227          77 :   N = gmul(ldata_get_conductor(ldata1), ldata_get_conductor(ldata2));
     228          77 :   Vga = shallowconcat(ldata_get_gammavec(ldata1), ldata_get_gammavec(ldata2));
     229          77 :   Vga = sort(Vga);
     230          77 :   eno = gmul(ldata_get_rootno(ldata1), ldata_get_rootno(ldata2));
     231          77 :   a1a2 = lfunconvol(ldata_get_an(ldata1), ldata_get_an(ldata2));
     232          77 :   b1b2 = lfuncombdual(lfunconvol, ldata1, ldata2);
     233          77 :   return mkvecn(r? 7: 6, a1a2, b1b2, Vga, k, N, eno, r);
     234             : }
     235             : 
     236             : GEN
     237          63 : lfunmul(GEN ldata1, GEN ldata2, long bitprec)
     238             : {
     239          63 :   pari_sp ltop = avma;
     240             :   GEN k;
     241          63 :   long prec = nbits2prec(bitprec);
     242          63 :   ldata1 = ldata_newprec(lfunmisc_to_ldata_shallow(ldata1), prec);
     243          63 :   ldata2 = ldata_newprec(lfunmisc_to_ldata_shallow(ldata2), prec);
     244          63 :   k = ldata_get_k(ldata1);
     245          63 :   if (!gequal(ldata_get_k(ldata2),k))
     246           0 :     pari_err_OP("lfunmul [weight]",ldata1, ldata2);
     247          63 :   return gerepilecopy(ltop, lfunmul_k(ldata1, ldata2, k, bitprec));
     248             : }
     249             : 
     250             : static GEN
     251          56 : lfundivpoles(GEN ldata1, GEN ldata2, long bitprec)
     252             : {
     253             :   long i, j, l;
     254          56 :   GEN k  = ldata_get_k(ldata1);
     255          56 :   GEN r1 = ldata_get_residue(ldata1);
     256          56 :   GEN r2 = ldata_get_residue(ldata2), r;
     257             : 
     258          56 :   if (r1 && typ(r1) != t_VEC) r1 = mkvec(mkvec2(k, r1));
     259          56 :   if (r2 && typ(r2) != t_VEC) r2 = mkvec(mkvec2(k, r2));
     260          56 :   if (!r1) return NULL;
     261          56 :   r1 = lfunrtopoles(r1);
     262          56 :   l = lg(r1); r = cgetg(l, t_VEC);
     263         112 :   for (i = j = 1; j < l; j++)
     264             :   {
     265          56 :     GEN be = gel(r1,j);
     266          56 :     GEN z = gdiv(lfun(ldata1,be,bitprec), lfun(ldata2,be,bitprec));
     267          56 :     if (valp(z) < 0) gel(r,i++) = mkvec2(be, z);
     268             :   }
     269          56 :   if (i == 1) return NULL;
     270          14 :   setlg(r, i); return r;
     271             : }
     272             : 
     273             : GEN
     274          56 : lfundiv(GEN ldata1, GEN ldata2, long bitprec)
     275             : {
     276          56 :   pari_sp ltop = avma;
     277             :   GEN k, r, N, v, v1, v2, eno, a1a2, b1b2, LD, eno2;
     278             :   long j, j1, j2, l1, l2;
     279          56 :   long prec = nbits2prec(bitprec);
     280          56 :   ldata1 = ldata_newprec(lfunmisc_to_ldata_shallow(ldata1), prec);
     281          56 :   ldata2 = ldata_newprec(lfunmisc_to_ldata_shallow(ldata2), prec);
     282          56 :   k = ldata_get_k(ldata1);
     283          56 :   if (!gequal(ldata_get_k(ldata2),k))
     284           0 :     pari_err_OP("lfundiv [weight]",ldata1, ldata2);
     285          56 :   r = lfundivpoles(ldata1, ldata2, bitprec);
     286          56 :   N = gdiv(ldata_get_conductor(ldata1), ldata_get_conductor(ldata2));
     287          56 :   if (typ(N) != t_INT) pari_err_OP("lfundiv [conductor]",ldata1, ldata2);
     288          56 :   a1a2 = lfunconvolinv(ldata_get_an(ldata1), ldata_get_an(ldata2));
     289          56 :   b1b2 = lfuncombdual(lfunconvolinv, ldata1, ldata2);
     290          56 :   eno2 = ldata_get_rootno(ldata2);
     291          56 :   eno = isintzero(eno2)? gen_0: gdiv(ldata_get_rootno(ldata1), eno2);
     292          56 :   v1 = shallowcopy(ldata_get_gammavec(ldata1));
     293          56 :   v2 = ldata_get_gammavec(ldata2);
     294          56 :   l1 = lg(v1); l2 = lg(v2);
     295         119 :   for (j2 = 1; j2 < l2; j2++)
     296             :   {
     297          91 :     for (j1 = 1; j1 < l1; j1++)
     298          91 :       if (gel(v1,j1) && gequal(gel(v1,j1), gel(v2,j2)))
     299             :       {
     300          63 :         gel(v1,j1) = NULL; break;
     301             :       }
     302          63 :     if (j1 == l1) pari_err_OP("lfundiv [Vga]",ldata1, ldata2);
     303             :   }
     304          56 :   v = cgetg(l1-l2+1, t_VEC);
     305         259 :   for (j1 = j = 1; j1 < l1; j1++)
     306         203 :     if (gel(v1, j1)) gel(v,j++) = gel(v1,j1);
     307             : 
     308          56 :   LD = mkvecn(7, a1a2, b1b2, v, k, N, eno, r);
     309          56 :   if (!r) setlg(LD,7);
     310          56 :   return gerepilecopy(ltop, LD);
     311             : }
     312             : 
     313             : static GEN
     314         245 : gamma_imagchi(GEN gam, GEN w)
     315             : {
     316         245 :   long i, j, k=1, l;
     317         245 :   GEN g = cgetg_copy(gam, &l);
     318         245 :   gam = shallowcopy(gam);
     319         735 :   for (i = l-1; i>=1; i--)
     320             :   {
     321         490 :     GEN al = gel(gam, i);
     322         490 :     if (al)
     323             :     {
     324         252 :       GEN N = gadd(w,gmul2n(real_i(al),1));
     325         252 :       if (gcmpgs(N,2) > 0)
     326             :       {
     327         238 :         GEN bl = gsubgs(al, 1);
     328         238 :         for (j=1; j < i; j++)
     329         238 :           if (gel(gam,j) && gequal(gel(gam,j), bl))
     330         238 :           { gel(gam,j) = NULL; break; }
     331         238 :         if (j==i) return NULL;
     332         238 :         gel(g, k++) = al;
     333         238 :         gel(g, k++) = bl;
     334          14 :       } else if (gequal0(N))
     335          14 :         gel(g, k++) = gaddgs(al, 1);
     336           0 :       else if (gequal1(N))
     337           0 :         gel(g, k++) = gsubgs(al, 1);
     338           0 :       else return NULL;
     339             :     }
     340             :   }
     341         245 :   return sort(g);
     342             : }
     343             : 
     344             : GEN
     345         889 : lfuntwist(GEN ldata1, GEN chi, long bitprec)
     346             : {
     347         889 :   pari_sp ltop = avma;
     348             :   GEN k, L, N, N1, N2, a, a1, a2, b, b1, b2, gam, gam1, gam2;
     349             :   GEN ldata2;
     350             :   long d1, t;
     351         889 :   long prec = nbits2prec(bitprec);
     352         889 :   ldata1 = ldata_newprec(lfunmisc_to_ldata_shallow(ldata1), prec);
     353         889 :   ldata2 = lfunmisc_to_ldata_shallow(chi);
     354         889 :   t = ldata_get_type(ldata2);
     355         889 :   a1 = ldata_get_an(ldata1);
     356         889 :   a2 = ldata_get_an(ldata2);
     357         889 :   if (t == t_LFUN_ZETA)
     358         427 :     return gerepilecopy(ltop, ldata1);
     359         462 :   if (t != t_LFUN_CHIZ && t != t_LFUN_KRONECKER &&
     360           7 :     ( t != t_LFUN_CHIGEN || nf_get_degree(bnr_get_nf(gmael(a2,2,1))) != 1))
     361           0 :     pari_err_TYPE("lfuntwist", chi);
     362         462 :   N1 = ldata_get_conductor(ldata1);
     363         462 :   N2 = ldata_get_conductor(ldata2);
     364         462 :   if (!gequal1(gcdii(N1, N2)))
     365           0 :     pari_err_IMPL("lfuntwist (conductors not coprime)");
     366         462 :   k = ldata_get_k(ldata1);
     367         462 :   d1 = ldata_get_degree(ldata1);
     368         462 :   N = gmul(N1, gpowgs(N2, d1));
     369         462 :   gam1 = ldata_get_gammavec(ldata1);
     370         462 :   gam2 = ldata_get_gammavec(ldata2);
     371         462 :   if (gequal0(gel(gam2, 1)))
     372         217 :     gam = gam1;
     373             :   else
     374         245 :     gam = gamma_imagchi(ldata_get_gammavec(ldata1), gaddgs(k,-1));
     375         462 :   if (!gam) pari_err_IMPL("lfuntwist (gammafactors)");
     376         462 :   b1 = ldata_get_dual(ldata1);
     377         462 :   b2 = ldata_get_dual(ldata2);
     378         462 :   a = tag(mkvec2(a1, a2), t_LFUN_TWIST);
     379         462 :   if (typ(b1)==t_INT)
     380         462 :     b = signe(b1) && signe(b2) ? gen_0: gen_1;
     381             :   else
     382           0 :     b = tag(mkvec2(b1,lfunconj(a2)), t_LFUN_TWIST);
     383         462 :   L = mkvecn(6, a, b, gam, k, N, gen_0);
     384         462 :   return gerepilecopy(ltop, L);
     385             : }
     386             : 
     387             : static GEN
     388         238 : lfundualpoles(GEN ldata, GEN reno)
     389             : {
     390             :   long l, j;
     391         238 :   GEN k = ldata_get_k(ldata);
     392         238 :   GEN r = gel(reno,2), eno = gel(reno,3), R;
     393         238 :   R = cgetg_copy(r, &l);
     394         728 :   for (j = 1; j < l; j++)
     395             :   {
     396         490 :     GEN b = gmael(r,j,1), e = gmael(r,j,2);
     397         490 :     long v = varn(e);
     398         490 :     GEN E = gsubst(gdiv(e, eno), v, gneg(pol_x(v)));
     399         490 :     gel(R,l-j) = mkvec2(gsub(k,b), E);
     400             :   }
     401         238 :   return R;
     402             : }
     403             : 
     404             : static GEN
     405         511 : ginvvec(GEN x)
     406             : {
     407         511 :   if (is_vec_t(typ(x)))
     408          42 :     pari_APPLY_same(ginv(gel(x,i)))
     409             :   else
     410         497 :     return ginv(x);
     411             : }
     412             : 
     413             : GEN
     414         560 : lfundual(GEN L, long bitprec)
     415             : {
     416         560 :   pari_sp av = avma;
     417         560 :   long prec = nbits2prec(bitprec);
     418         560 :   GEN ldata = ldata_newprec(lfunmisc_to_ldata_shallow(L), prec);
     419         560 :   GEN a = ldata_get_an(ldata), b = ldata_get_dual(ldata);
     420         560 :   GEN e = ldata_get_rootno(ldata);
     421         560 :   GEN ldual, ad, bd, ed, Rd = NULL;
     422         560 :   if (typ(b) == t_INT)
     423             :   {
     424         532 :     ad = equali1(b) ? lfunconj(a): a;
     425         532 :     bd = b;
     426             :   }
     427          28 :   else { ad = b; bd = a; }
     428         560 :   if (lg(ldata)==8)
     429             :   {
     430         238 :     GEN reno = lfunrootres(ldata, bitprec);
     431         238 :     e = gel(reno,3);
     432         238 :     Rd = lfundualpoles(ldata, reno);
     433             :   }
     434         560 :   ed = isintzero(e) ? e: ginvvec(e);
     435         560 :   ldual = mkvecn(Rd ? 7:6, ad, bd, gel(ldata,3), gel(ldata,4), gel(ldata,5), ed, Rd);
     436         560 :   return gerepilecopy(av, ldual);
     437             : }
     438             : 
     439             : static GEN
     440          49 : RgV_Rg_translate(GEN x, GEN s)
     441         140 : { pari_APPLY_same(gadd(gel(x,i),s)) }
     442             : 
     443             : static GEN
     444          28 : pole_translate(GEN x, GEN s, GEN Ns)
     445             : {
     446          28 :   x = shallowcopy(x);
     447          28 :   gel(x,1) = gadd(gel(x,1), s);
     448          28 :   if (Ns)
     449          28 :     gel(x,2) = gmul(gel(x,2), Ns);
     450          28 :   return x;
     451             : }
     452             : 
     453             : static GEN
     454          14 : poles_translate(GEN x, GEN s, GEN Ns)
     455          42 : { pari_APPLY_same(pole_translate(gel(x,i), s, Ns)) }
     456             : 
     457             : /* r / x + O(1) */
     458             : static GEN
     459         217 : simple_pole(GEN r)
     460             : {
     461             :   GEN S;
     462         217 :   if (isintzero(r)) return gen_0;
     463         189 :   S = deg1ser_shallow(gen_0, r, 1);
     464         189 :   setvalp(S, -1); return S;
     465             : }
     466             : 
     467             : GEN
     468          49 : lfunshift(GEN ldata, GEN s, long flag, long bitprec)
     469             : {
     470          49 :   pari_sp ltop = avma;
     471             :   GEN k, k1, L, N, a, b, gam, eps, res;
     472          49 :   long prec = nbits2prec(bitprec);
     473          49 :   if (!is_rational_t(typ(s))) pari_err_TYPE("lfunshift",s);
     474          49 :   ldata = ldata_newprec(lfunmisc_to_ldata_shallow(ldata), prec);
     475          49 :   a = ldata_get_an(ldata);
     476          49 :   b = ldata_get_dual(ldata);
     477          49 :   gam = RgV_Rg_translate(ldata_get_gammavec(ldata), gneg(s));
     478          49 :   k = gadd(ldata_get_k(ldata), gmul2n(s, 1));
     479          49 :   k1 = gadd(ldata_get_k1(ldata), s);
     480          49 :   N = ldata_get_conductor(ldata);
     481          49 :   eps = ldata_get_rootno(ldata);
     482          49 :   res = ldata_get_residue(ldata);
     483          49 :   a = tag(mkvec2(a, s), t_LFUN_SHIFT);
     484          49 :   if (typ(b) != t_INT)
     485           0 :     b = tag(mkvec2(b, s), t_LFUN_SHIFT);
     486          49 :   if (res)
     487          49 :     switch(typ(res))
     488             :     {
     489           0 :     case t_VEC:
     490           0 :       res = poles_translate(res, s, NULL);
     491           0 :       break;
     492          14 :     case t_COL:
     493          14 :       res = poles_translate(res, s, gpow(N, gmul2n(s, -1), prec));
     494          14 :       break;
     495          35 :     default:
     496          35 :       res = mkvec(mkvec2(gsub(k, s), simple_pole(res)));
     497             :     }
     498          49 :   L = mkvecn(res ? 7: 6, a, b, gam, mkvec2(k, k1), N, eps, res);
     499          49 :   if (flag) L = lfunmul_k(ldata, L, gsub(k, s), bitprec);
     500          49 :   return gerepilecopy(ltop, L);
     501             : }
     502             : 
     503             : /*****************************************************************/
     504             : /*  L-series from closure                                        */
     505             : /*****************************************************************/
     506             : static GEN
     507       56805 : localfactor(void *E, GEN p, long n)
     508             : {
     509       56805 :   GEN s = closure_callgen2((GEN)E, p, utoi(n));
     510       56805 :   return direuler_factor(s, n);
     511             : }
     512             : static GEN
     513        1890 : vecan_closure(GEN a, long L, long prec)
     514             : {
     515        1890 :   long ta = typ(a);
     516        1890 :   GEN gL, Sbad = NULL;
     517             : 
     518        1890 :   if (!L) return cgetg(1,t_VEC);
     519        1890 :   if (ta == t_VEC)
     520             :   {
     521         973 :     long l = lg(a);
     522         973 :     if (l == 1) pari_err_TYPE("vecan_closure", a);
     523         973 :     ta = typ(gel(a,1));
     524             :     /* regular vector, return it */
     525         973 :     if (ta != t_CLOSURE) return vecslice(a, 1, minss(L,l-1));
     526         119 :     if (l != 3) pari_err_TYPE("vecan_closure", a);
     527         119 :     Sbad = gel(a,2);
     528         119 :     if (typ(Sbad) != t_VEC) pari_err_TYPE("vecan_closure", a);
     529         112 :     a = gel(a,1);
     530             :   }
     531         917 :   else if (ta != t_CLOSURE) pari_err_TYPE("vecan_closure", a);
     532        1015 :   push_localprec(prec);
     533        1015 :   gL = stoi(L);
     534        1015 :   switch(closure_arity(a))
     535             :   {
     536         371 :     case 2:
     537         371 :       a = direuler_bad((void*)a, localfactor, gen_2, gL,gL, Sbad);
     538         336 :       break;
     539         637 :     case 1:
     540         637 :       a = closure_callgen1(a, gL);
     541         637 :       if (typ(a) != t_VEC) pari_err_TYPE("vecan_closure", a);
     542         630 :       break;
     543           7 :     default: pari_err_TYPE("vecan_closure [wrong arity]", a);
     544             :       a = NULL; /*LCOV_EXCL_LINE*/
     545             :   }
     546         966 :   pop_localprec(); return a;
     547             : }
     548             : 
     549             : /*****************************************************************/
     550             : /*  L-series of Dirichlet characters.                            */
     551             : /*****************************************************************/
     552             : 
     553             : static GEN
     554        1652 : lfunzeta(void)
     555             : {
     556        1652 :   GEN zet = mkvecn(7, NULL, gen_0, NULL, gen_1, gen_1, gen_1, gen_1);
     557        1652 :   gel(zet,1) = tag(gen_1, t_LFUN_ZETA);
     558        1652 :   gel(zet,3) = mkvec(gen_0);
     559        1652 :   return zet;
     560             : }
     561             : static GEN
     562         182 : lfunzetainit(GEN dom, long der, long bitprec)
     563         182 : { return lfuninit(lfunzeta(), dom, der, bitprec); }
     564             : 
     565             : static GEN
     566        1288 : vecan_Kronecker(GEN D, long n)
     567             : {
     568        1288 :   GEN v = cgetg(n+1, t_VECSMALL);
     569        1288 :   ulong Du = itou_or_0(D);
     570        1288 :   long i, id, d = Du ? minuu(Du, n): n;
     571       31850 :   for (i = 1; i <= d; i++) v[i] = krois(D,i);
     572      162197 :   for (id = i; i <= n; i++,id++) /* periodic mod d */
     573             :   {
     574      160909 :     if (id > d) id = 1;
     575      160909 :     gel(v, i) = gel(v, id);
     576             :   }
     577        1288 :   return v;
     578             : }
     579             : 
     580             : static GEN
     581        4018 : lfunchiquad(GEN D)
     582             : {
     583             :   GEN r;
     584        4018 :   D = coredisc(D);
     585        4018 :   if (equali1(D)) return lfunzeta();
     586        3598 :   if (!isfundamental(D)) pari_err_TYPE("lfunchiquad [not primitive]", D);
     587        3598 :   r = mkvecn(6, NULL, gen_0, NULL, gen_1, NULL, gen_1);
     588        3598 :   gel(r,1) = tag(icopy(D), t_LFUN_KRONECKER);
     589        3598 :   gel(r,3) = mkvec(signe(D) < 0? gen_1: gen_0);
     590        3598 :   gel(r,5) = mpabs(D);
     591        3598 :   return r;
     592             : }
     593             : 
     594             : /* Begin Hecke characters. Here a character is assumed to be given by a
     595             :    vector on the generators of the ray class group clgp of CL_m(K).
     596             :    If clgp = [h,[d1,...,dk],[g1,...,gk]] with dk|...|d2|d1, a character chi
     597             :    is given by [a1,a2,...,ak] such that chi(gi)=\zeta_di^ai. */
     598             : 
     599             : /* Value of CHI on x, coprime to bnr.mod */
     600             : static GEN
     601       81179 : chigeneval_i(GEN logx, GEN d, GEN nchi, GEN z, long prec)
     602             : {
     603       81179 :   pari_sp av = avma;
     604       81179 :   GEN e = FpV_dotproduct(nchi, logx, d);
     605       81179 :   if (!is_vec_t(typ(z)))
     606        1120 :     return gerepileupto(av, gpow(z, e, prec));
     607             :   else
     608             :   {
     609       80059 :     ulong i = itou(e);
     610       80059 :     set_avma(av); return gel(z, i+1);
     611             :   }
     612             : }
     613             : 
     614             : static GEN
     615       76895 : chigenevalvec(GEN logx, GEN nchi, GEN z, long prec, long multi)
     616             : {
     617       76895 :   GEN d = gel(nchi,1), x = gel(nchi, 2);
     618       76895 :   if (multi)
     619       12852 :     pari_APPLY_same(chigeneval_i(logx, d, gel(x,i), z, prec))
     620             :   else
     621       72611 :     return chigeneval_i(logx, d, x, z, prec);
     622             : }
     623             : 
     624             : /* return x + yz; y != 0; z = 0,1 "often"; x = 0 "often" */
     625             : static GEN
     626     1781598 : gaddmul(GEN x, GEN y, GEN z)
     627             : {
     628             :   pari_sp av;
     629     1781598 :   if (typ(z) == t_INT)
     630             :   {
     631     1595104 :     if (!signe(z)) return x;
     632       26005 :     if (equali1(z)) return gadd(x,y);
     633             :   }
     634      205366 :   if (isintzero(x)) return gmul(y,z);
     635      113253 :   av = avma;
     636      113253 :   return gerepileupto(av, gadd(x, gmul(y,z)));
     637             : }
     638             : 
     639             : static GEN
     640     1750574 : gaddmulvec(GEN x, GEN y, GEN z, long multi)
     641             : {
     642     1750574 :   if (multi)
     643       93072 :     pari_APPLY_same(gaddmul(gel(x,i),gel(y,i),gel(z,i)))
     644             :   else
     645     1719550 :     return gaddmul(x,y,z);
     646             : }
     647             : 
     648             : static GEN
     649        1897 : mkvchi(GEN chi, long n)
     650             : {
     651             :   GEN v;
     652        1897 :   if (lg(chi) > 1 && is_vec_t(typ(gel(chi,1))))
     653         245 :   {
     654         245 :     long d = lg(chi)-1;
     655         245 :     v = const_vec(n, zerovec(d));
     656         245 :     gel(v,1) = const_vec(d, gen_1);
     657             :   }
     658             :   else
     659        1652 :     v = vec_ei(n, 1);
     660        1897 :   return v;
     661             : }
     662             : 
     663             : static GEN
     664         987 : vecan_chiZ(GEN an, long n, long prec)
     665             : {
     666             :   forprime_t iter;
     667         987 :   GEN G = gel(an,1);
     668         987 :   GEN nchi = gel(an,2), gord = gel(nchi,1), chi = gel(nchi,2), z;
     669         987 :   GEN gp = cgetipos(3), v = mkvchi(chi, n);
     670         987 :   GEN N = znstar_get_N(G);
     671         987 :   long ord = itos_or_0(gord);
     672         987 :   ulong Nu = itou_or_0(N);
     673         987 :   long i, id, d = Nu ? minuu(Nu, n): n;
     674         987 :   long multichi= (lg(chi) > 1 && is_vec_t(typ(gel(chi,1))));
     675             :   ulong p;
     676         987 :   if (!multichi && ord && n > (ord>>4))
     677         861 :   {
     678         861 :     GEN w = ncharvecexpo(G, nchi);
     679         861 :     z = grootsof1(ord, prec);
     680       14777 :     for (i = 1; i <= d; i++)
     681       13916 :       if (w[i] >= 0) gel(v, i) = gel(z, w[i]+1);
     682             :   }
     683             :   else
     684             :   {
     685         126 :     z = rootsof1_cx(gord, prec);
     686         126 :     u_forprime_init(&iter, 2, d);
     687         805 :     while ((p = u_forprime_next(&iter)))
     688             :     {
     689             :       GEN ch;
     690             :       ulong k;
     691         679 :       if (!umodiu(N,p)) continue;
     692         560 :       gp[2] = p;
     693         560 :       ch = chigenevalvec(znconreylog(G, gp), nchi, z, prec, multichi);
     694         560 :       gel(v, p)  = ch;
     695        1582 :       for (k = 2*p; k <= (ulong)d; k += p)
     696        1022 :         gel(v, k) = gaddmulvec(gel(v, k), ch, gel(v, k/p), multichi);
     697             :     }
     698             :   }
     699      254639 :   for (id = i = d+1; i <= n; i++,id++) /* periodic mod d */
     700             :   {
     701      253652 :     if (id > d) id = 1;
     702      253652 :     gel(v, i) = gel(v, id);
     703             :   }
     704         987 :   return v;
     705             : }
     706             : 
     707             : static GEN
     708         910 : vecan_chigen(GEN an, long n, long prec)
     709             : {
     710             :   forprime_t iter;
     711         910 :   GEN bnr = gel(an,1), nf = bnr_get_nf(bnr);
     712         910 :   GEN nchi = gel(an,2), gord = gel(nchi,1), chi = gel(nchi,2), z;
     713         910 :   GEN gp = cgetipos(3), v = mkvchi(chi, n);
     714         910 :   GEN N = gel(bnr_get_mod(bnr), 1), NZ = gcoeff(N,1,1);
     715         910 :   long ord = itos_or_0(gord);
     716         910 :   long multichi= (lg(chi) > 1 && is_vec_t(typ(gel(chi,1))));
     717             :   ulong p;
     718             : 
     719         910 :   if (ord && n > (ord>>4))
     720         910 :     z = grootsof1(ord, prec);
     721             :   else
     722           0 :     z = rootsof1_cx(gord, prec);
     723             : 
     724         910 :   if (nf_get_degree(nf) == 1)
     725             :   {
     726         616 :     ulong Nu = itou_or_0(NZ);
     727         616 :     long i, id, d = Nu ? minuu(Nu, n): n;
     728         616 :     u_forprime_init(&iter, 2, d);
     729        3066 :     while ((p = u_forprime_next(&iter)))
     730             :     {
     731             :       GEN ch;
     732             :       ulong k;
     733        2450 :       if (!umodiu(NZ,p)) continue;
     734        1771 :       gp[2] = p;
     735        1771 :       ch = chigenevalvec(isprincipalray(bnr,gp), nchi, z, prec, multichi);
     736        1771 :       gel(v, p)  = ch;
     737        3906 :       for (k = 2*p; k <= (ulong)d; k += p)
     738        2135 :         gel(v, k) = gaddmulvec(gel(v, k), ch, gel(v, k/p), multichi);
     739             :     }
     740        8414 :     for (id = i = d+1; i <= n; i++,id++) /* periodic mod d */
     741             :     {
     742        7798 :       if (id > d) id = 1;
     743        7798 :       gel(v, i) = gel(v, id);
     744             :     }
     745             :   }
     746             :   else
     747             :   {
     748         294 :     GEN BOUND = stoi(n);
     749         294 :     u_forprime_init(&iter, 2, n);
     750       75061 :     while ((p = u_forprime_next(&iter)))
     751             :     {
     752             :       GEN L;
     753             :       long j;
     754       74767 :       int check = !umodiu(NZ,p);
     755       74767 :       gp[2] = p;
     756       74767 :       L = idealprimedec_limit_norm(nf, gp, BOUND);
     757      149450 :       for (j = 1; j < lg(L); j++)
     758             :       {
     759       74683 :         GEN pr = gel(L, j), ch;
     760             :         ulong k, q;
     761       74683 :         if (check && idealval(nf, N, pr)) continue;
     762       74564 :         ch = chigenevalvec(isprincipalray(bnr,pr), nchi, z, prec, multichi);
     763       74564 :         q = upr_norm(pr);
     764       74564 :         gel(v, q) = gadd(gel(v, q), ch);
     765     1821981 :         for (k = 2*q; k <= (ulong)n; k += q)
     766     1747417 :           gel(v, k) = gaddmulvec(gel(v, k), ch, gel(v, k/q), multichi);
     767             :       }
     768             :     }
     769             :   }
     770         910 :   return v;
     771             : }
     772             : 
     773             : static GEN lfunzetak_i(GEN T);
     774             : static GEN
     775        3941 : vec01(long r1, long r2)
     776             : {
     777        3941 :   long d = r1+r2, i;
     778        3941 :   GEN v = cgetg(d+1,t_VEC);
     779       10493 :   for (i = 1; i <= r1; i++) gel(v,i) = gen_0;
     780        6083 :   for (     ; i <= d;  i++) gel(v,i) = gen_1;
     781        3941 :   return v;
     782             : }
     783             : /* v = vector of normalized characters of order dividing o; renormalize
     784             :  * so that all have same apparent order o */
     785             : static GEN
     786          21 : char_renormalize(GEN v, GEN o)
     787             : {
     788             :   long i, l;
     789          21 :   GEN w = cgetg_copy(v, &l);
     790          63 :   for (i = 1; i < l; i++)
     791             :   {
     792          42 :     GEN C = gel(v,i), oc = gel(C,1), c = gel(C,2);
     793          42 :     if (!equalii(o, oc)) c = gmul(c, diviiexact(o, oc));
     794          42 :     gel(w,i) = c;
     795             :   }
     796          21 :   return w;
     797             : }
     798             : /* G is a bid of nftyp typ_BIDZ */
     799             : static GEN
     800        1722 : lfunchiZ(GEN G, GEN CHI)
     801             : {
     802        1722 :   pari_sp av = avma;
     803        1722 :   GEN sig = NULL, N = bid_get_ideal(G), nchi, r;
     804             :   int real;
     805             :   long s;
     806             : 
     807        1722 :   if (typ(N) != t_INT) pari_err_TYPE("lfunchiZ", G);
     808        1722 :   if (typ(CHI) == t_VEC && !RgV_is_ZV(CHI))
     809          14 :   {
     810          35 :     GEN C, G0 = G, o = gen_1;
     811          35 :     long i, l = lg(CHI);
     812          35 :     nchi = cgetg(l, t_VEC);
     813          35 :     N = znconreyconductor(G, gel(CHI,1), &C);
     814          28 :     if (typ(N) != t_INT) G = znstar0(N, 1);
     815          28 :     s = zncharisodd(G, C);
     816          70 :     for (i = 1; i < l; i++)
     817             :     {
     818          56 :       if (i > 1)
     819             :       {
     820          28 :         if (!gequal(N, znconreyconductor(G0, gel(CHI,i), &C))
     821          21 :             || zncharisodd(G, C) != s)
     822          14 :           pari_err_TYPE("lfuncreate [different conductors]", CHI);
     823             :       }
     824          42 :       C = znconreylog_normalize(G, C);
     825          42 :       o = lcmii(o, gel(C,1)); /* lcm with charorder */
     826          42 :       gel(nchi,i) = C;
     827             :     }
     828          14 :     nchi = mkvec2(o, char_renormalize(nchi, o));
     829          14 :     if (typ(N) != t_INT) N = gel(N,1);
     830             :   }
     831             :   else
     832             :   {
     833        1687 :     N = znconreyconductor(G, CHI, &CHI);
     834        1687 :     if (typ(N) != t_INT)
     835             :     {
     836           0 :       if (equali1(gel(N,1))) { set_avma(av); return lfunzeta(); }
     837           0 :       G = znstar0(N, 1);
     838           0 :       N = gel(N,1);
     839             :     }
     840             :     /* CHI now primitive on G */
     841        1687 :     switch(itou_or_0(zncharorder(G, CHI)))
     842             :     {
     843         427 :       case 1: set_avma(av); return lfunzeta();
     844         658 :       case 2: if (zncharisodd(G,CHI)) N = negi(N);
     845         658 :               return gerepileupto(av, lfunchiquad(N));
     846             :     }
     847         602 :     nchi = znconreylog_normalize(G, CHI);
     848         602 :     s = zncharisodd(G, CHI);
     849             :   }
     850         616 :   sig = mkvec(s? gen_1: gen_0);
     851         616 :   real = abscmpiu(gel(nchi,1), 2) <= 0;
     852         616 :   r = mkvecn(6, tag(mkvec2(G,nchi), t_LFUN_CHIZ),
     853             :                 real? gen_0: gen_1, sig, gen_1, N, gen_0);
     854         616 :   return gerepilecopy(av, r);
     855             : }
     856             : 
     857             : static GEN
     858        1064 : lfunchigen(GEN bnr, GEN CHI)
     859             : {
     860        1064 :   pari_sp av = avma;
     861             :   GEN N, sig, Ldchi, nf, nchi, NN;
     862             :   long r1, r2, n1;
     863             :   int real;
     864             : 
     865        1064 :   if (typ(CHI) == t_VEC && !RgV_is_ZV(CHI))
     866           7 :   {
     867          14 :     long map, i, l = lg(CHI);
     868          14 :     GEN chi = gel(CHI,1), bnr0 = bnr, o = gen_1;
     869          14 :     GEN D = cyc_normalize(bnr_get_cyc(bnr));
     870          14 :     nchi = cgetg(l, t_VEC);
     871          14 :     N = bnr_get_mod(bnr);
     872          14 :     bnr_char_sanitize(&bnr, &chi);
     873          14 :     map = (bnr != bnr0);
     874          35 :     for (i = 1; i < l; i++)
     875             :     {
     876          28 :       if (i > 1)
     877             :       {
     878          14 :         chi = gel(CHI,i);
     879          14 :         if (!map)
     880             :         {
     881          14 :           if (!bnrisconductor(bnr, chi))
     882           7 :             pari_err_TYPE("lfuncreate [different conductors]", CHI);
     883             :         }
     884             :         else
     885             :         {
     886           0 :           if (!gequal(bnrconductor_raw(bnr0, chi), N))
     887           0 :             pari_err_TYPE("lfuncreate [different conductors]", CHI);
     888           0 :           chi = bnrchar_primitive_raw(bnr0, bnr, chi);
     889             :         }
     890             :       }
     891          21 :       chi = char_normalize(chi, D);
     892          21 :       o = lcmii(o, gel(chi,1)); /* lcm with charorder */
     893          21 :       gel(nchi,i) = chi;
     894             :     }
     895           7 :     nchi = mkvec2(o, char_renormalize(nchi, o));
     896             :   }
     897             :   else
     898             :   {
     899        1050 :     bnr_char_sanitize(&bnr, &CHI);
     900        1050 :     nchi = NULL; /* now CHI is primitive wrt bnr */
     901             :   }
     902             : 
     903        1057 :   N = bnr_get_mod(bnr);
     904        1057 :   nf = bnr_get_nf(bnr);
     905        1057 :   n1 = lg(vec01_to_indices(gel(N,2))) - 1; /* vecsum(N[2]) */
     906        1057 :   N = gel(N,1);
     907        1057 :   NN = mulii(idealnorm(nf, N), absi_shallow(nf_get_disc(nf)));
     908        1057 :   if (!nchi)
     909             :   {
     910        1050 :     if (equali1(NN)) { set_avma(av); return lfunzeta(); }
     911         644 :     if (ZV_equal0(CHI)) return gerepilecopy(av, lfunzetak_i(bnr_get_bnf(bnr)));
     912         637 :     nchi = char_normalize(CHI, cyc_normalize(bnr_get_cyc(bnr)));
     913             :   }
     914         644 :   real = abscmpiu(gel(nchi,1), 2) <= 0;
     915         644 :   nf_get_sign(nf, &r1, &r2);
     916         644 :   sig = vec01(r1+r2-n1, r2+n1);
     917         644 :   Ldchi = mkvecn(6, tag(mkvec2(bnr, nchi), t_LFUN_CHIGEN),
     918             :                     real? gen_0: gen_1, sig, gen_1, NN, gen_0);
     919         644 :   return gerepilecopy(av, Ldchi);
     920             : }
     921             : 
     922             : /* Find all characters of clgp whose kernel contain group given by HNF H.
     923             :  * Set *pcnj[i] if chi[i] is not real */
     924             : static GEN
     925         413 : chigenkerfind(GEN bnr, GEN H, GEN *pcnj)
     926             : {
     927         413 :   GEN res, cnj, L = bnrchar(bnr, H, NULL), cyc = bnr_get_cyc(bnr);
     928         413 :   long i, k, l = lg(L);
     929             : 
     930         413 :   res = cgetg(l, t_VEC);
     931         413 :   *pcnj = cnj = cgetg(l, t_VECSMALL);
     932        1519 :   for (i = k = 1; i < l; i++)
     933             :   {
     934        1106 :     GEN chi = gel(L,i), c = charconj(cyc, chi);
     935        1106 :     long fl = ZV_cmp(c, chi);
     936        1106 :     if (fl < 0) continue; /* keep one char in pair of conjugates */
     937         931 :     gel(res, k) = chi;
     938         931 :     cnj[k] = fl; k++;
     939             :   }
     940         413 :   setlg(cnj, k);
     941         413 :   setlg(res, k); return res;
     942             : }
     943             : 
     944             : static GEN
     945        1400 : lfunzetak_i(GEN T)
     946             : {
     947        1400 :   GEN Vga, N, nf, bnf = checkbnf_i(T), r = gen_0/*unknown*/;
     948             :   long r1, r2;
     949             : 
     950        1400 :   if (bnf)
     951         105 :     nf = bnf_get_nf(bnf);
     952             :   else
     953             :   {
     954        1295 :     nf = checknf_i(T);
     955        1295 :     if (!nf) nf = T = nfinit(T, DEFAULTPREC);
     956             :   }
     957        1400 :   nf_get_sign(nf,&r1,&r2);
     958        1400 :   N = absi_shallow(nf_get_disc(nf));
     959        1400 :   if (bnf)
     960             :   {
     961         105 :     GEN h = bnf_get_no(bnf);
     962         105 :     GEN R = bnf_get_reg(bnf);
     963         105 :     long prec = nf_get_prec(nf);
     964         210 :     r = gdiv(gmul(mulir(shifti(h, r1+r2), powru(mppi(prec), r2)), R),
     965         105 :              mulur(bnf_get_tuN(bnf), gsqrt(N, prec)));
     966             :   }
     967        1400 :   Vga = vec01(r1+r2,r2);
     968        1400 :   return mkvecn(7, tag(T,t_LFUN_NF), gen_0, Vga, gen_1, N, gen_1, r);
     969             : }
     970             : static GEN
     971         616 : lfunzetak(GEN T)
     972         616 : { pari_sp ltop = avma; return gerepilecopy(ltop, lfunzetak_i(T)); }
     973             : 
     974             : /* bnf = NULL: base field = Q */
     975             : GEN
     976         413 : lfunabelianrelinit(GEN nfabs, GEN bnf, GEN polrel, GEN dom, long der, long bitprec)
     977             : {
     978         413 :   pari_sp ltop = avma;
     979             :   GEN cond, chi, cnj, res, bnr, M, domain;
     980             :   long l, i;
     981         413 :   long v = -1;
     982             : 
     983         413 :   if (bnf) bnf = checkbnf(bnf);
     984             :   else
     985             :   {
     986         406 :     v = fetch_var();
     987         406 :     bnf = Buchall(pol_x(v), 0, nbits2prec(bitprec));
     988             :   }
     989         413 :   if (typ(polrel) != t_POL) pari_err_TYPE("lfunabelianrelinit", polrel);
     990         413 :   cond = rnfconductor0(bnf, polrel, 1);
     991         413 :   bnr = gel(cond,2);
     992         413 :   chi = chigenkerfind(bnr, gel(cond,3), &cnj);
     993         413 :   l = lg(chi); res = cgetg(l, t_VEC);
     994        1344 :   for (i = 1; i < l; ++i)
     995             :   {
     996         931 :     GEN L = lfunchigen(bnr, gel(chi,i));
     997         931 :     gel(res, i) = lfuninit(L, dom, der, bitprec);
     998             :   }
     999         413 :   if (v >= 0) delete_var();
    1000         413 :   M = mkvec3(res, const_vecsmall(l-1, 1), cnj);
    1001         413 :   domain = mkvec2(dom, mkvecsmall2(der, bitprec));
    1002         413 :   return gerepilecopy(ltop, lfuninit_make(t_LDESC_PRODUCT, lfunzetak_i(nfabs), M, domain));
    1003             : }
    1004             : 
    1005             : /*****************************************************************/
    1006             : /*                 Dedekind zeta functions                       */
    1007             : /*****************************************************************/
    1008             : static GEN
    1009        2163 : dirzetak0(GEN nf, ulong N)
    1010             : {
    1011        2163 :   GEN vect, c, c2, T = nf_get_pol(nf), index = nf_get_index(nf);
    1012        2163 :   pari_sp av = avma, av2;
    1013        2163 :   const ulong SQRTN = usqrt(N);
    1014             :   ulong i, p, lx;
    1015        2163 :   long court[] = {evaltyp(t_INT)|_evallg(3), evalsigne(1)|evallgefint(3),0};
    1016             :   forprime_t S;
    1017             : 
    1018        2163 :   c  = cgetalloc(t_VECSMALL, N+1);
    1019        2163 :   c2 = cgetalloc(t_VECSMALL, N+1);
    1020     3272122 :   c2[1] = c[1] = 1; for (i=2; i<=N; i++) c[i] = 0;
    1021        2163 :   u_forprime_init(&S, 2, N); av2 = avma;
    1022      419167 :   while ( (p = u_forprime_next(&S)) )
    1023             :   {
    1024      417004 :     set_avma(av2);
    1025      417004 :     if (umodiu(index, p)) /* p does not divide index */
    1026      416724 :       vect = gel(Flx_degfact(ZX_to_Flx(T,p), p),1);
    1027             :     else
    1028             :     {
    1029         280 :       court[2] = p;
    1030         280 :       vect = idealprimedec_degrees(nf,court);
    1031             :     }
    1032      417004 :     lx = lg(vect);
    1033      417004 :     if (p <= SQRTN)
    1034       37464 :       for (i=1; i<lx; i++)
    1035             :       {
    1036       25459 :         ulong qn, q = upowuu(p, vect[i]); /* Norm P[i] */
    1037       25459 :         if (!q || q > N) break;
    1038       21763 :         memcpy(c2 + 2, c + 2, (N-1)*sizeof(long));
    1039             :         /* c2[i] <- c[i] + sum_{k = 1}^{v_q(i)} c[i/q^k] for all i <= N */
    1040       44849 :         for (qn = q; qn <= N; qn *= q)
    1041             :         {
    1042       44849 :           ulong k0 = N/qn, k, k2; /* k2 = k*qn */
    1043     4064102 :           for (k = k0, k2 = k*qn; k > 0; k--, k2 -=qn) c2[k2] += c[k];
    1044       44849 :           if (q > k0) break; /* <=> q*qn > N */
    1045             :         }
    1046       21763 :         swap(c, c2);
    1047             :       }
    1048             :     else /* p > sqrt(N): simpler */
    1049      794444 :       for (i=1; i<lx; i++)
    1050             :       {
    1051             :         ulong k, k2; /* k2 = k*p */
    1052      704004 :         if (vect[i] > 1) break;
    1053             :         /* c2[i] <- c[i] + sum_{k = 1}^{v_q(i)} c[i/q^k] for all i <= N */
    1054     2346722 :         for (k = N/p, k2 = k*p; k > 0; k--, k2 -= p) c[k2] += c[k];
    1055             :       }
    1056             :   }
    1057        2163 :   set_avma(av);
    1058        2163 :   pari_free(c2); return c;
    1059             : }
    1060             : 
    1061             : GEN
    1062        2163 : dirzetak(GEN nf, GEN b)
    1063             : {
    1064             :   GEN z, c;
    1065             :   long n;
    1066             : 
    1067        2163 :   if (typ(b) != t_INT) pari_err_TYPE("dirzetak",b);
    1068        2163 :   if (signe(b) <= 0) return cgetg(1,t_VEC);
    1069        2163 :   nf = checknf(nf);
    1070        2163 :   n = itou_or_0(b); if (!n) pari_err_OVERFLOW("dirzetak");
    1071        2163 :   c = dirzetak0(nf, n);
    1072        2163 :   z = vecsmall_to_vec(c); pari_free(c); return z;
    1073             : }
    1074             : 
    1075             : static GEN
    1076         686 : linit_get_mat(GEN linit)
    1077             : {
    1078         686 :   if (linit_get_type(linit)==t_LDESC_PRODUCT)
    1079         161 :     return lfunprod_get_fact(linit_get_tech(linit));
    1080             :   else
    1081         525 :     return mkvec3(mkvec(linit), mkvecsmall(1), mkvecsmall(0));
    1082             : }
    1083             : 
    1084             : static GEN
    1085         343 : lfunproduct(GEN ldata, GEN linit1, GEN linit2, GEN domain)
    1086             : {
    1087         343 :   GEN M1 = linit_get_mat(linit1);
    1088         343 :   GEN M2 = linit_get_mat(linit2);
    1089         343 :   GEN M3 = mkvec3(shallowconcat(gel(M1, 1), gel(M2, 1)),
    1090         343 :                   vecsmall_concat(gel(M1, 2), gel(M2, 2)),
    1091         343 :                   vecsmall_concat(gel(M1, 3), gel(M2, 3)));
    1092         343 :   return lfuninit_make(t_LDESC_PRODUCT, ldata, M3, domain);
    1093             : }
    1094             : 
    1095             : static GEN
    1096           0 : lfunzetakinit_raw(GEN T, GEN dom, long der, long bitprec)
    1097             : {
    1098           0 :   pari_sp ltop = avma;
    1099           0 :   GEN ldata = lfunzetak_i(T);
    1100           0 :   return gerepileupto(ltop, lfuninit(ldata, dom, der, bitprec));
    1101             : }
    1102             : 
    1103             : static GEN
    1104         343 : lfunzetakinit_quotient(GEN nf, GEN polk, GEN dom, long der, long bitprec)
    1105             : {
    1106         343 :   pari_sp av = avma;
    1107             :   GEN ak, an, nfk, Vga, ldata, N, Lk, LKk, domain;
    1108             :   long r1k, r2k, r1, r2;
    1109             : 
    1110         343 :   nf_get_sign(nf,&r1,&r2);
    1111         343 :   nfk = nfinit(polk, nbits2prec(bitprec));
    1112         343 :   Lk = lfunzetakinit(nfk, dom, der, 0, bitprec); /* zeta_k */
    1113         343 :   nf_get_sign(nfk,&r1k,&r2k);
    1114         343 :   Vga = vec01((r1+r2) - (r1k+r2k), r2-r2k);
    1115         343 :   N = absi_shallow(diviiexact(nf_get_disc(nf), nf_get_disc(nfk)));
    1116         343 :   ak = nf_get_degree(nf)==1 ? tag(gen_1, t_LFUN_ZETA): tag(nfk, t_LFUN_NF);
    1117         343 :   an = tag(mkvec2(tag(nf,t_LFUN_NF), ak), t_LFUN_DIV);
    1118         343 :   ldata = mkvecn(6, an, gen_0, Vga, gen_1, N, gen_1);
    1119         343 :   LKk = lfuninit(ldata, dom, der, bitprec); /* zeta_K/zeta_k */
    1120         343 :   domain = mkvec2(dom, mkvecsmall2(der, bitprec));
    1121         343 :   return gerepilecopy(av, lfunproduct(lfunzetak_i(nf), Lk, LKk, domain));
    1122             : }
    1123             : 
    1124             : static GEN
    1125             : lfunzetakinit_artin(GEN nf, GEN gal, GEN dom, long der, long bitprec);
    1126             : 
    1127             : static GEN
    1128         427 : lfunzetakinit_Galois(GEN nf, GEN gal, GEN dom, long der, long bitprec)
    1129             : {
    1130         427 :   GEN grp = galois_group(gal);
    1131         427 :   if (group_isabelian(grp))
    1132         406 :     return lfunabelianrelinit(nf, NULL, gal_get_pol(gal), dom, der, bitprec);
    1133          21 :   else return lfunzetakinit_artin(nf, gal, dom, der, bitprec);
    1134             : }
    1135             : 
    1136             : GEN
    1137         952 : lfunzetakinit(GEN NF, GEN dom, long der, long flag, long bitprec)
    1138             : {
    1139         952 :   GEN nf = checknf(NF);
    1140             :   GEN G, nfs, sbg;
    1141         952 :   long lf, d = nf_get_degree(nf);
    1142         952 :   if (d == 1) return lfunzetainit(dom, der, bitprec);
    1143         770 :   G = galoisinit(nf, NULL);
    1144         770 :   if (!isintzero(G))
    1145         427 :     return lfunzetakinit_Galois(nf, G, dom, der, bitprec);
    1146         343 :   nfs = nfsubfields(nf, 0); lf = lg(nfs)-1;
    1147         343 :   sbg = gmael(nfs,lf-1,1);
    1148         343 :   if (flag && d > 4*degpol(sbg))
    1149           0 :     return lfunzetakinit_raw(nf, dom, der, bitprec);
    1150         343 :   return lfunzetakinit_quotient(nf, sbg, dom, der, bitprec);
    1151             : }
    1152             : 
    1153             : /***************************************************************/
    1154             : /*             Elliptic Curves and Modular Forms               */
    1155             : /***************************************************************/
    1156             : 
    1157             : static GEN
    1158         175 : lfunellnf(GEN e)
    1159             : {
    1160         175 :   pari_sp av = avma;
    1161         175 :   GEN ldata = cgetg(7, t_VEC), nf = ellnf_get_nf(e);
    1162         175 :   GEN N = gel(ellglobalred(e), 1);
    1163         175 :   long n = nf_get_degree(nf);
    1164         175 :   gel(ldata, 1) = tag(e, t_LFUN_ELL);
    1165         175 :   gel(ldata, 2) = gen_0;
    1166         175 :   gel(ldata, 3) = vec01(n, n);
    1167         175 :   gel(ldata, 4) = gen_2;
    1168         175 :   gel(ldata, 5) = mulii(idealnorm(nf,N), sqri(nf_get_disc(nf)));
    1169         175 :   gel(ldata, 6) = stoi(ellrootno_global(e));
    1170         175 :   return gerepilecopy(av, ldata);
    1171             : }
    1172             : 
    1173             : static GEN
    1174        1561 : lfunellQ(GEN e)
    1175             : {
    1176        1561 :   pari_sp av = avma;
    1177        1561 :   GEN ldata = cgetg(7, t_VEC);
    1178        1561 :   gel(ldata, 1) = tag(ellanal_globalred(e, NULL), t_LFUN_ELL);
    1179        1561 :   gel(ldata, 2) = gen_0;
    1180        1561 :   gel(ldata, 3) = mkvec2(gen_0, gen_1);
    1181        1561 :   gel(ldata, 4) = gen_2;
    1182        1561 :   gel(ldata, 5) = ellQ_get_N(e);
    1183        1561 :   gel(ldata, 6) = stoi(ellrootno_global(e));
    1184        1561 :   return gerepilecopy(av, ldata); /* ellanal_globalred not gerepile-safe */
    1185             : }
    1186             : 
    1187             : static GEN
    1188        1736 : lfunell(GEN e)
    1189             : {
    1190        1736 :   long t = ell_get_type(e);
    1191        1736 :   switch(t)
    1192             :   {
    1193        1561 :     case t_ELL_Q: return lfunellQ(e);
    1194         175 :     case t_ELL_NF:return lfunellnf(e);
    1195             :   }
    1196           0 :   pari_err_TYPE("lfun",e);
    1197             :   return NULL; /*LCOV_EXCL_LINE*/
    1198             : }
    1199             : 
    1200             : static GEN
    1201         140 : ellsympow_gamma(long m)
    1202             : {
    1203         140 :   GEN V = cgetg(m+2, t_VEC);
    1204         140 :   long i = 1, j;
    1205         140 :   if (!odd(m)) gel(V, i++) = stoi(-2*(m>>2));
    1206         364 :   for (j = (m+1)>>1; j > 0; i+=2, j--)
    1207             :   {
    1208         224 :     gel(V,i)   = stoi(1-j);
    1209         224 :     gel(V,i+1) = stoi(1-j+1);
    1210             :   }
    1211         140 :   return V;
    1212             : }
    1213             : 
    1214             : static GEN
    1215       78475 : ellsympow_trace(GEN p, GEN t, long m)
    1216             : {
    1217       78475 :   long k, n = m >> 1;
    1218       78475 :   GEN tp = gpowers0(sqri(t), n, odd(m)? t: NULL);
    1219       78495 :   GEN pp = gen_1, b = gen_1, r = gel(tp,n+1);
    1220      215583 :   for(k=1; k<=n; k++)
    1221             :   {
    1222             :     GEN s;
    1223      137341 :     pp = mulii(pp, p);
    1224      136409 :     b  = diviuexact(muliu(b, (m-(2*k-1))*(m-(2*k-2))), k*(m-(k-1)));
    1225      136244 :     s = mulii(mulii(b, gel(tp,1+n-k)), pp);
    1226      136956 :     r = odd(k) ? subii(r, s): addii(r, s);
    1227             :   }
    1228       78242 :   return r;
    1229             : }
    1230             : 
    1231             : static GEN
    1232        3059 : ellsympow_abelian(GEN p, GEN ap, long m, long o)
    1233             : {
    1234        3059 :   pari_sp av = avma;
    1235        3059 :   long i, M, n = (m+1)>>1;
    1236             :   GEN pk, tv, pn, pm, F, v;
    1237        3059 :   if (!odd(o))
    1238             :   {
    1239           0 :     if (odd(m)) return pol_1(0);
    1240           0 :     M = m >> 1; o >>= 1;
    1241             :   }
    1242             :   else
    1243        3059 :     M = m * ((o+1) >> 1);
    1244        3059 :   pk = gpowers(p,n); pn = gel(pk,n+1);
    1245        3059 :   tv = cgetg(m+2,t_VEC);
    1246        3059 :   gel(tv, 1) = gen_2;
    1247        3059 :   gel(tv, 2) = ap;
    1248       10318 :   for (i = 3; i <= m+1; i++)
    1249        7259 :     gel(tv,i) = subii(mulii(ap,gel(tv,i-1)), mulii(p,gel(tv,i-2)));
    1250        3059 :   pm = odd(m)? mulii(gel(pk,n), pn): sqri(pn); /* cheap p^m */
    1251        3059 :   F = deg2pol_shallow(pm, gen_0, gen_1, 0);
    1252        3059 :   v = odd(m) ? pol_1(0): deg1pol_shallow(negi(pn), gen_1, 0);
    1253        9177 :   for (i = M % o; i < n; i += o) /* o | m-2*i */
    1254             :   {
    1255        6118 :     gel(F,3) = negi(mulii(gel(tv,m-2*i+1), gel(pk,i+1)));
    1256        6118 :     v = ZX_mul(v, F);
    1257             :   }
    1258        3059 :   return gerepilecopy(av, v);
    1259             : }
    1260             : 
    1261             : static GEN
    1262       81565 : ellsympow(GEN E, ulong m, GEN p, long n)
    1263             : {
    1264       81565 :   pari_sp av = avma;
    1265       81565 :   GEN ap = ellap(E, p);
    1266       81519 :   if (n <= 2)
    1267             :   {
    1268       78463 :     GEN t = ellsympow_trace(p, ap, m);
    1269       78255 :     return deg1pol_shallow(t, gen_1, 0);
    1270             :   }
    1271             :   else
    1272        3056 :     return gerepileupto(av, RgXn_inv_i(ellsympow_abelian(p, ap, m, 1), n));
    1273             : }
    1274             : 
    1275             : GEN
    1276        5467 : direllsympow_worker(GEN P, ulong X, GEN E, ulong m)
    1277             : {
    1278        5467 :   pari_sp av = avma;
    1279        5467 :   long i, l = lg(P);
    1280        5467 :   GEN W = cgetg(l, t_VEC);
    1281       87030 :   for(i = 1; i < l; i++)
    1282             :   {
    1283       81563 :     ulong p = uel(P,i);
    1284       81563 :     long d = ulogint(X, p) + 1; /* minimal d such that p^d > X */
    1285       81577 :     gel(W,i) = ellsympow(E, m, utoi(uel(P,i)), d);
    1286             :   }
    1287        5467 :   return gerepilecopy(av, mkvec2(P,W));
    1288             : }
    1289             : 
    1290             : static GEN
    1291         343 : vecan_ellsympow(GEN an, long n)
    1292             : {
    1293         343 :   GEN nn = utoi(n), crvm = gel(an,1), bad = gel(an,2);
    1294         343 :   GEN worker = snm_closure(is_entry("_direllsympow_worker"), crvm);
    1295         343 :   return pardireuler(worker, gen_2, nn, nn, bad);
    1296             : }
    1297             : 
    1298             : static long
    1299         196 : ellsympow_betam(long o, long m)
    1300             : {
    1301         196 :   const long c3[]={3, -1, 1};
    1302         196 :   const long c12[]={6, -2, 2, 0, 4, -4};
    1303         196 :   const long c24[]={12, -2, -4, 6, 4, -10};
    1304         196 :   if (!odd(o) && odd(m)) return 0;
    1305         161 :   switch(o)
    1306             :   {
    1307           0 :     case 1:  return m+1;
    1308          14 :     case 2:  return m+1;
    1309          84 :     case 3:  case 6: return (m+c3[m%3])/3;
    1310           0 :     case 4:  return m%4 == 0 ? (m+2)/2: m/2;
    1311          21 :     case 8:  return m%4 == 0 ? (m+4)/4: (m-2)/4;
    1312          35 :     case 12: return (m+c12[(m%12)/2])/6;
    1313           7 :     case 24: return (m+c24[(m%12)/2])/12;
    1314             :   }
    1315           0 :   return 0;
    1316             : }
    1317             : 
    1318             : static long
    1319          98 : ellsympow_epsm(long o, long m) { return m + 1 - ellsympow_betam(o, m); }
    1320             : 
    1321             : static GEN
    1322          98 : ellsympow_multred(GEN E, GEN p, long m, long vN, long *cnd, long *w)
    1323             : {
    1324          98 :   if (vN == 1 || !odd(m))
    1325             :   {
    1326          98 :     GEN s = (odd(m) && signe(ellap(E,p)) < 0)? gen_1: gen_m1;
    1327          98 :     *cnd = m;
    1328          98 :     *w = odd(m)? ellrootno(E, p): 1;
    1329          98 :     return deg1pol_shallow(s, gen_1, 0);
    1330             :   }
    1331             :   else
    1332             :   {
    1333           0 :     *cnd = equaliu(p,2)? ((m+1)>>1) * vN: m+1;
    1334           0 :     *w = (m & 3) == 1? ellrootno(E, p): 1;
    1335           0 :     return pol_1(0);
    1336             :   }
    1337             : }
    1338             : 
    1339             : static GEN
    1340          98 : ellsympow_nonabelian(GEN p, long m, long bet)
    1341             : {
    1342          98 :  GEN q = powiu(p, m >> 1), q2 = sqri(q), F;
    1343          98 :  if (odd(m))
    1344             :  {
    1345          35 :    q2 = mulii(q2, p); /* p^m */
    1346          35 :    return gpowgs(deg2pol_shallow(q2, gen_0, gen_1, 0), bet>>1);
    1347             :  }
    1348          63 :  togglesign_safe(&q2);
    1349          63 :  F = gpowgs(deg2pol_shallow(q2, gen_0, gen_1, 0), bet>>1);
    1350          63 :  if (!odd(bet)) return F;
    1351          28 :  if (m%4 != 2) togglesign_safe(&q);
    1352          28 :  return gmul(F, deg1pol_shallow(q, gen_1, 0));
    1353             : }
    1354             : 
    1355             : static long
    1356           0 : safe_Z_pvalrem(GEN n, GEN p, GEN *pr)
    1357           0 : { return signe(n)==0? -1: Z_pvalrem(n, p, pr); }
    1358             : 
    1359             : static GEN
    1360           0 : c4c6_ap(GEN c4, GEN c6, GEN p)
    1361             : {
    1362           0 :   GEN N = Fp_ellcard(Fp_muls(c4, -27, p), Fp_muls(c6, -54, p), p);
    1363           0 :   return subii(addiu(p, 1), N);
    1364             : }
    1365             : 
    1366             : static GEN
    1367           0 : ellsympow_abelian_twist(GEN E, GEN p, long m, long o)
    1368             : {
    1369           0 :   GEN ap, c4t, c6t, c4 = ell_get_c4(E), c6 = ell_get_c6(E);
    1370           0 :   long v4 = safe_Z_pvalrem(c4, p, &c4t);
    1371           0 :   long v6 = safe_Z_pvalrem(c6, p, &c6t);
    1372           0 :   if (v6>=0 && (v4==-1 || 3*v4>=2*v6)) c6 = c6t;
    1373           0 :   if (v4>=0 && (v6==-1 || 3*v4<=2*v6)) c4 = c4t;
    1374           0 :   ap = c4c6_ap(c4, c6, p);
    1375           0 :   return ellsympow_abelian(p, ap, m, o);
    1376             : }
    1377             : 
    1378             : static GEN
    1379           0 : ellsympow_goodred(GEN E, GEN p, long m, long *cnd, long *w)
    1380             : {
    1381           0 :   long o = 12/cgcd(12, Z_pval(ell_get_disc(E), p));
    1382           0 :   long bet = ellsympow_betam(o, m);
    1383           0 :   long eps = m + 1 - bet;
    1384           0 :   *w = odd(m) && odd(eps>>1) ? ellrootno(E,p): 1;
    1385           0 :   *cnd = eps;
    1386           0 :   if (umodiu(p, o) == 1)
    1387           0 :     return ellsympow_abelian_twist(E, p, m, o);
    1388             :   else
    1389           0 :     return ellsympow_nonabelian(p, m, bet);
    1390             : }
    1391             : 
    1392             : static long
    1393          70 : ellsympow_inertia3(GEN E, long vN)
    1394             : {
    1395          70 :   long vD = Z_lval(ell_get_disc(E), 3);
    1396          70 :   if (vN==2) return vD%2==0 ? 2: 4;
    1397          70 :   if (vN==4) return vD%4==0 ? 3: 6;
    1398          70 :   if (vN==3 || vN==5) return 12;
    1399           0 :   return 0;
    1400             : }
    1401             : 
    1402             : static long
    1403          70 : ellsympow_deltam3(long o, long m, long vN)
    1404             : {
    1405          70 :   if (o==3 || o==6) return ellsympow_epsm(3, m);
    1406          70 :   if (o==12 && vN ==3) return (ellsympow_epsm(3, m))/2;
    1407           0 :   if (o==12 && vN ==5) return (ellsympow_epsm(3, m))*3/2;
    1408           0 :   return 0;
    1409             : }
    1410             : 
    1411             : static long
    1412           0 : ellsympow_isabelian3(GEN E)
    1413             : {
    1414           0 :   ulong c4 = umodiu(ell_get_c4(E),81), c6 = umodiu(ell_get_c6(E), 243);
    1415           0 :   return (c4 == 27 || (c4%27==9 && (c6==108 || c6==135)));
    1416             : }
    1417             : 
    1418             : static long
    1419          35 : ellsympow_rootno3(GEN E, GEN p, long o, long m)
    1420             : {
    1421          35 :   const long  w6p[]={1,-1,-1,-1,1,1};
    1422          35 :   const long  w6n[]={-1,1,-1,1,-1,1};
    1423          35 :   const long w12p[]={1,1,-1,1,1,1};
    1424          35 :   const long w12n[]={-1,-1,-1,-1,-1,1};
    1425          35 :   long w = ellrootno(E, p), mm = (m%12)>>1;
    1426          35 :   switch(o)
    1427             :   {
    1428           0 :     case 2: return m%4== 1 ? -1: 1;
    1429           0 :     case 6:  return w == 1 ? w6p[mm]: w6n[mm];
    1430          35 :     case 12: return w == 1 ? w12p[mm]: w12n[mm];
    1431           0 :     default: return 1;
    1432             :   }
    1433             : }
    1434             : 
    1435             : static GEN
    1436          70 : ellsympow_goodred3(GEN E, GEN F, GEN p, long m, long vN, long *cnd, long *w)
    1437             : {
    1438          70 :   long o = ellsympow_inertia3(E, vN);
    1439          70 :   long bet = ellsympow_betam(o, m);
    1440          70 :   *cnd = m + 1 - bet + ellsympow_deltam3(o, m, vN);
    1441          70 :   *w = odd(m)? ellsympow_rootno3(E, p, o, m): 1;
    1442          70 :   if (o==1 || o==2)
    1443           0 :     return ellsympow_abelian(p, ellap(F, p), m, o);
    1444          70 :   if ((o==3 || o==6) && ellsympow_isabelian3(F))
    1445           0 :     return ellsympow_abelian(p, p, m, o);
    1446             :   else
    1447          70 :     return ellsympow_nonabelian(p, m, bet);
    1448             : }
    1449             : 
    1450             : static long
    1451          28 : ellsympow_inertia2(GEN F, long vN)
    1452             : {
    1453          28 :   long vM = itos(gel(elllocalred(F, gen_2),1));
    1454          28 :   GEN c6 = ell_get_c6(F);
    1455          28 :   long v6 = signe(c6) ? vali(c6): 24;
    1456          28 :   if (vM==0) return vN==0 ? 1: 2;
    1457          28 :   if (vM==2) return vN==2 ? 3: 6;
    1458          14 :   if (vM==5) return 8;
    1459           7 :   if (vM==8) return v6>=9? 8: 4;
    1460           7 :   if (vM==3 || vN==7) return 24;
    1461           0 :   return 0;
    1462             : }
    1463             : 
    1464             : static long
    1465          28 : ellsympow_deltam2(long o, long m, long vN)
    1466             : {
    1467          28 :   if ((o==2 || o==6) && vN==4) return ellsympow_epsm(2, m);
    1468          28 :   if ((o==2 || o==6) && vN==6) return 2*ellsympow_epsm(2, m);
    1469          28 :   if (o==4) return 2*ellsympow_epsm(4, m)+ellsympow_epsm(2, m);
    1470          28 :   if (o==8 && vN==5) return ellsympow_epsm(8, m)+ellsympow_epsm(2, m)/2;
    1471          21 :   if (o==8 && vN==6) return ellsympow_epsm(8, m)+ellsympow_epsm(2, m);
    1472          21 :   if (o==8 && vN==8) return ellsympow_epsm(8, m)+ellsympow_epsm(4, m)+ellsympow_epsm(2, m);
    1473          21 :   if (o==24 && vN==3) return (2*ellsympow_epsm(8, m)+ellsympow_epsm(2, m))/6;
    1474          14 :   if (o==24 && vN==4) return (ellsympow_epsm(8, m)+ellsympow_epsm(2, m)*2)/3;
    1475          14 :   if (o==24 && vN==6) return (ellsympow_epsm(8, m)+ellsympow_epsm(2, m)*5)/3;
    1476          14 :   if (o==24 && vN==7) return (ellsympow_epsm(8, m)*10+ellsympow_epsm(2, m)*5)/6;
    1477          14 :   return 0;
    1478             : }
    1479             : 
    1480             : static long
    1481           0 : ellsympow_isabelian2(GEN F)
    1482           0 : { return umodi2n(ell_get_c4(F),7) == 96; }
    1483             : 
    1484             : static long
    1485           0 : ellsympow_rootno2(GEN E, long vN, long m, long bet)
    1486             : {
    1487           0 :   long eps2 = (m + 1 - bet)>>1;
    1488           0 :   long eta = odd(vN) && m%8==3 ? -1 : 1;
    1489           0 :   long w2 = odd(eps2) ? ellrootno(E, gen_2): 1;
    1490           0 :   return eta == w2 ? 1 : -1;
    1491             : }
    1492             : 
    1493             : static GEN
    1494          28 : ellsympow_goodred2(GEN E, GEN F, GEN p, long m, long vN, long *cnd, long *w)
    1495             : {
    1496          28 :   long o = ellsympow_inertia2(F, vN);
    1497          28 :   long bet = ellsympow_betam(o, m);
    1498          28 :   *cnd = m + 1 - bet + ellsympow_deltam2(o, m, vN);
    1499          28 :   *w = odd(m) ? ellsympow_rootno2(E, vN, m, bet): 1;
    1500          28 :   if (o==1 || o==2)
    1501           0 :     return ellsympow_abelian(p, ellap(F, p), m, o);
    1502          28 :   if (o==4 && ellsympow_isabelian2(F))
    1503           0 :     return ellsympow_abelian(p, p, m, o);
    1504             :   else
    1505          28 :     return ellsympow_nonabelian(p, m, bet);
    1506             : }
    1507             : 
    1508             : static GEN
    1509         189 : ellminimaldotwist(GEN E, GEN *pD, long prec)
    1510             : {
    1511         189 :   GEN D = ellminimaltwistcond(E), Et = elltwist(E, D), Etmin;
    1512         189 :   if (pD) *pD = D;
    1513         189 :   Et = ellinit(Et, NULL, prec);
    1514         189 :   Etmin = ellminimalmodel(Et, NULL);
    1515         189 :   obj_free(Et); return Etmin;
    1516             : }
    1517             : 
    1518             : /* Based on
    1519             : Symmetric powers of elliptic curve L-functions,
    1520             : Phil Martin and Mark Watkins, ANTS VII
    1521             : <http://magma.maths.usyd.edu.au/users/watkins/papers/antsVII.pdf>
    1522             : with thanks to Mark Watkins. BA20180402
    1523             : */
    1524             : static GEN
    1525         140 : lfunellsympow(GEN e, ulong m)
    1526             : {
    1527         140 :   pari_sp av = avma;
    1528             :   GEN B, N, Nfa, pr, ex, ld, bad, ejd, et, pole;
    1529         140 :   long i, l, mero, w = (m&7)==1 || (m&7)==3 ? -1: 1;
    1530         140 :   checkell_Q(e);
    1531         140 :   e = ellminimalmodel(e, NULL);
    1532         140 :   ejd = Q_denom(ell_get_j(e));
    1533         140 :   mero = m==0 || (m%4==0 && ellQ_get_CM(e)<0);
    1534         140 :   ellQ_get_Nfa(e, &N, &Nfa);
    1535         140 :   pr = gel(Nfa,1);
    1536         140 :   ex = gel(Nfa,2); l = lg(pr);
    1537         140 :   if (ugcd(umodiu(N,6), 6) == 1)
    1538          21 :     et = NULL;
    1539             :   else
    1540         119 :     et = ellminimaldotwist(e, NULL, DEFAULTPREC);
    1541         140 :   B = gen_1;
    1542         140 :   bad = cgetg(l, t_VEC);
    1543         336 :   for (i=1; i<l; i++)
    1544             :   {
    1545         196 :     long vN = itos(gel(ex,i));
    1546         196 :     GEN p = gel(pr,i), eul;
    1547             :     long cnd, wp;
    1548         196 :     if (dvdii(ejd, p))
    1549          98 :       eul = ellsympow_multred(e, p, m, vN, &cnd, &wp);
    1550          98 :     else if (equaliu(p, 2))
    1551          28 :       eul = ellsympow_goodred2(e, et, p, m, vN, &cnd, &wp);
    1552          70 :     else if (equaliu(p, 3))
    1553          70 :       eul = ellsympow_goodred3(e, et, p, m, vN, &cnd, &wp);
    1554             :     else
    1555           0 :       eul = ellsympow_goodred(e, p, m, &cnd, &wp);
    1556         196 :     gel(bad, i) = mkvec2(p, ginv(eul));
    1557         196 :     B = mulii(B, powiu(p,cnd));
    1558         196 :     w *= wp;
    1559             :   }
    1560         140 :   pole = mero ? mkvec(mkvec2(stoi(1+(m>>1)),gen_0)): NULL;
    1561         280 :   ld = mkvecn(mero? 7: 6, tag(mkvec2(mkvec2(e,utoi(m)),bad), t_LFUN_SYMPOW_ELL),
    1562         140 :         gen_0, ellsympow_gamma(m), stoi(m+1), B, stoi(w), pole);
    1563         140 :   if (et) obj_free(et);
    1564         140 :   return gerepilecopy(av, ld);
    1565             : }
    1566             : 
    1567             : GEN
    1568          70 : lfunsympow(GEN ldata, ulong m)
    1569             : {
    1570          70 :   ldata = lfunmisc_to_ldata_shallow(ldata);
    1571          70 :   if (ldata_get_type(ldata) != t_LFUN_ELL)
    1572           0 :     pari_err_IMPL("lfunsympow");
    1573          70 :   return lfunellsympow(gel(ldata_get_an(ldata), 2), m);
    1574             : }
    1575             : 
    1576             : static GEN
    1577          28 : lfunmfspec_i(GEN lmisc, long bit)
    1578             : {
    1579             :   GEN linit, ldataf, v, ve, vo, om, op, B, dom;
    1580             :   long k, k2, j;
    1581             : 
    1582          28 :   ldataf = lfunmisc_to_ldata_shallow(lmisc);
    1583          28 :   if (!gequal(ldata_get_gammavec(ldataf), mkvec2(gen_0,gen_1)))
    1584           0 :     pari_err_TYPE("lfunmfspec", lmisc);
    1585          28 :   k = gtos(ldata_get_k(ldataf));
    1586          28 :   if (k == 1) return mkvec2(cgetg(1, t_VEC), gen_1);
    1587          21 :   dom = mkvec3(dbltor(k/2.), dbltor((k-2)/2.), gen_0);
    1588          21 :   if (is_linit(lmisc) && linit_get_type(lmisc) == t_LDESC_INIT
    1589           0 :       && sdomain_isincl((double)k, dom, lfun_get_dom(linit_get_tech(lmisc))))
    1590           0 :     linit = lmisc;
    1591             :   else
    1592          21 :     linit = lfuninit(ldataf, dom, 0, bit);
    1593          21 :   B = int2n(bit/4);
    1594          21 :   v = cgetg(k, t_VEC);
    1595         168 :   for (j = 1; j < k; j++) gel(v,j) = lfunlambda(linit, utoi(j), bit);
    1596          21 :   om = gel(v,1);
    1597          21 :   if (odd(k)) return mkvec2(bestappr(gdiv(v, om), B), om);
    1598             : 
    1599           7 :   k2 = k/2;
    1600           7 :   ve = cgetg(k2, t_VEC);
    1601           7 :   vo = cgetg(k2+1, t_VEC);
    1602           7 :   gel(vo,1) = om;
    1603          42 :   for (j = 1; j < k2; j++)
    1604             :   {
    1605          35 :     gel(ve,j) = gel(v,2*j);
    1606          35 :     gel(vo,j+1) = gel(v,2*j+1);
    1607             :   }
    1608           7 :   if (k2 == 1) { om = gen_1;    op = gel(v,1); }
    1609           7 :   else         { om = gel(v,2); op = gel(v,3); }
    1610           7 :   if (maxss(gexpo(imag_i(om)), gexpo(imag_i(op))) > -bit/2)
    1611           0 :     pari_err_TYPE("lfunmfspec", lmisc);
    1612           7 :   ve = gdiv(ve, om);
    1613           7 :   vo = gdiv(vo, op);
    1614           7 :   return mkvec4(bestappr(ve,B), bestappr(vo,B), om, op);
    1615             : }
    1616             : GEN
    1617          28 : lfunmfspec(GEN lmisc, long bit)
    1618             : {
    1619          28 :   pari_sp av = avma;
    1620          28 :   return gerepilecopy(av, lfunmfspec_i(lmisc, bit));
    1621             : }
    1622             : 
    1623             : static long
    1624          28 : ellsymsq_bad2(GEN c4, GEN c6, long e)
    1625             : {
    1626          28 :   switch (e)
    1627             :   {
    1628          14 :     case 2: return 1;
    1629           7 :     case 3: return 0;
    1630           7 :     case 5: return 0;
    1631           0 :     case 7: return 0;
    1632           0 :     case 8:
    1633           0 :       if (!umodi2n(c6,9)) return 0;
    1634           0 :       return umodi2n(c4,7)==32 ? 1 : -1;
    1635           0 :     default: return 0;
    1636             :   }
    1637             : }
    1638             : static long
    1639          14 : ellsymsq_bad3(GEN c4, GEN c6, long e)
    1640             : {
    1641             :   long c6_243, c4_81;
    1642          14 :   switch (e)
    1643             :   {
    1644           0 :     case 2: return 1;
    1645          14 :     case 3: return 0;
    1646           0 :     case 5: return 0;
    1647           0 :     case 4:
    1648           0 :       c4_81 = umodiu(c4,81);
    1649           0 :       if (c4_81 == 27) return -1;
    1650           0 :       if (c4_81%27 != 9) return 1;
    1651           0 :       c6_243 = umodiu(c6,243);
    1652           0 :       return (c6_243==108 || c6_243==135)? -1: 1;
    1653           0 :     default: return 0;
    1654             :   }
    1655             : }
    1656             : static int
    1657           0 : c4c6_testp(GEN c4, GEN c6, GEN p)
    1658           0 : { GEN p2 = sqri(p); return (dvdii(c6,p2) && !dvdii(c4,p2)); }
    1659             : /* assume e = v_p(N) >= 2 */
    1660             : static long
    1661          42 : ellsymsq_badp(GEN c4, GEN c6, GEN p, long e)
    1662             : {
    1663          42 :   if (absequaliu(p, 2)) return ellsymsq_bad2(c4, c6, e);
    1664          14 :   if (absequaliu(p, 3)) return ellsymsq_bad3(c4, c6, e);
    1665           0 :   switch(umodiu(p, 12UL))
    1666             :   {
    1667           0 :     case 1: return -1;
    1668           0 :     case 5: return c4c6_testp(c4,c6,p)? -1: 1;
    1669           0 :     case 7: return c4c6_testp(c4,c6,p)?  1:-1;
    1670           0 :     default:return 1; /* p%12 = 11 */
    1671             :   }
    1672             : }
    1673             : static GEN
    1674          70 : lfunellsymsqmintwist(GEN e)
    1675             : {
    1676          70 :   pari_sp av = avma;
    1677             :   GEN N, Nfa, P, E, V, c4, c6, ld;
    1678             :   long i, l, k;
    1679          70 :   checkell_Q(e);
    1680          70 :   e = ellminimalmodel(e, NULL);
    1681          70 :   ellQ_get_Nfa(e, &N, &Nfa);
    1682          70 :   c4 = ell_get_c4(e);
    1683          70 :   c6 = ell_get_c6(e);
    1684          70 :   P = gel(Nfa,1); l = lg(P);
    1685          70 :   E = gel(Nfa,2);
    1686          70 :   V = cgetg(l, t_VEC);
    1687         196 :   for (i=k=1; i<l; i++)
    1688             :   {
    1689         126 :     GEN p = gel(P,i);
    1690         126 :     long a, e = itos(gel(E,i));
    1691         126 :     if (e == 1) continue;
    1692          42 :     a = ellsymsq_badp(c4, c6, p, e);
    1693          42 :     gel(V,k++) = mkvec2(p, stoi(a));
    1694             :   }
    1695          70 :   setlg(V, k);
    1696          70 :   ld = lfunellsympow(e, 2);
    1697          70 :   return gerepilecopy(av, mkvec2(ld, V));
    1698             : }
    1699             : 
    1700             : static GEN
    1701          70 : mfpeters(GEN ldata2, GEN fudge, GEN N, long k, long bitprec)
    1702             : {
    1703          70 :   GEN t, L = real_i(lfun(ldata2, stoi(k), bitprec));
    1704          70 :   long prec = nbits2prec(bitprec);
    1705          70 :   t = powrs(mppi(prec), k+1); shiftr_inplace(t, 2*k-1); /* Pi/2 * (4Pi)^k */
    1706          70 :   return gmul(gdiv(gmul(mulii(N,mpfact(k-1)), fudge), t), L);
    1707             : }
    1708             : 
    1709             : /* Assume E to be twist-minimal */
    1710             : static GEN
    1711          70 : lfunellmfpetersmintwist(GEN E, long bitprec)
    1712             : {
    1713          70 :   pari_sp av = avma;
    1714          70 :   GEN symsq, veceuler, N = ellQ_get_N(E), fudge = gen_1;
    1715          70 :   long j, k = 2;
    1716          70 :   symsq = lfunellsymsqmintwist(E);
    1717          70 :   veceuler = gel(symsq,2);
    1718         112 :   for (j = 1; j < lg(veceuler); j++)
    1719             :   {
    1720          42 :     GEN v = gel(veceuler,j), p = gel(v,1), q = powis(p,1-k);
    1721          42 :     long s = signe(gel(v,2));
    1722          42 :     if (s) fudge = gmul(fudge, s==1 ? gaddsg(1, q): gsubsg(1, q));
    1723             :   }
    1724          70 :   return gerepileupto(av, mfpeters(gel(symsq,1),fudge,N,k,bitprec));
    1725             : }
    1726             : 
    1727             : /* From Christophe Delaunay, http://delaunay.perso.math.cnrs.fr/these.pdf */
    1728             : static GEN
    1729          70 : elldiscfix(GEN E, GEN Et, GEN D)
    1730             : {
    1731          70 :   GEN N = ellQ_get_N(E), Nt = ellQ_get_N(Et);
    1732          70 :   GEN P = gel(absZ_factor(D), 1);
    1733          70 :   GEN f = gen_1;
    1734          70 :   long i, l = lg(P);
    1735         133 :   for (i=1; i < l; i++)
    1736             :   {
    1737          63 :     GEN r, p = gel(P,i);
    1738          63 :     long v = Z_pval(N, p), vt = Z_pval(Nt, p);
    1739          63 :     if (v <= vt) continue;
    1740             :     /* v > vt */
    1741          49 :     if (absequaliu(p, 2))
    1742             :     {
    1743          28 :       if (vt == 0 && v >= 4)
    1744           0 :         r = shifti(subsi(9, sqri(ellap(Et, p))), v-3);  /* 9=(2+1)^2 */
    1745          28 :       else if (vt == 1)
    1746           7 :         r = gmul2n(utoipos(3), v-3);  /* not in Z if v=2 */
    1747          21 :       else if (vt >= 2)
    1748          21 :         r = int2n(v-vt);
    1749             :       else
    1750           0 :         r = gen_1; /* vt = 0, 1 <= v <= 3 */
    1751             :     }
    1752          21 :     else if (vt >= 1)
    1753          14 :       r = gdiv(subiu(sqri(p), 1), p);
    1754             :     else
    1755           7 :       r = gdiv(mulii(subiu(p, 1), subii(sqri(addiu(p, 1)), sqri(ellap(Et, p)))), p);
    1756          49 :     f = gmul(f, r);
    1757             :   }
    1758          70 :   return f;
    1759             : }
    1760             : 
    1761             : GEN
    1762          70 : lfunellmfpeters(GEN E, long bitprec)
    1763             : {
    1764          70 :   pari_sp ltop = avma;
    1765          70 :   GEN D, Et = ellminimaldotwist(E, &D, nbits2prec(bitprec));
    1766          70 :   GEN nor = lfunellmfpetersmintwist(Et, bitprec);
    1767          70 :   GEN nor2 = gmul(nor, elldiscfix(E, Et, D));
    1768          70 :   obj_free(Et); return gerepileupto(ltop, nor2);
    1769             : }
    1770             : 
    1771             : /*************************************************************/
    1772             : /*               Genus 2 curves                              */
    1773             : /*************************************************************/
    1774             : 
    1775             : static void
    1776      185919 : Flv_diffnext(GEN d, ulong p)
    1777             : {
    1778      185919 :   long j, n = lg(d)-1;
    1779     1298855 :   for(j = n; j>=2; j--)
    1780     1113119 :     uel(d,j) = Fl_add(uel(d,j), uel(d,j-1), p);
    1781      185736 : }
    1782             : 
    1783             : static GEN
    1784        1764 : Flx_difftable(GEN P, ulong p)
    1785             : {
    1786        1764 :   long i, n = degpol(P);
    1787        1764 :   GEN V = cgetg(n+2, t_VECSMALL);
    1788        1764 :   uel(V, n+1) = Flx_constant(P);
    1789       12345 :   for(i = n; i >= 1; i--)
    1790             :   {
    1791       10581 :     P = Flx_diff1(P, p);
    1792       10580 :     uel(V, i) = Flx_constant(P);
    1793             :   }
    1794        1764 :   return V;
    1795             : }
    1796             : 
    1797             : static long
    1798        1763 : Flx_genus2trace_naive(GEN H, ulong p)
    1799             : {
    1800        1763 :   pari_sp av = avma;
    1801             :   ulong i, j;
    1802        1763 :   long a, n = degpol(H);
    1803        1763 :   GEN k = const_vecsmall(p, -1), d;
    1804        1764 :   k[1] = 0;
    1805       96255 :   for (i=1, j=1; i < p; i += 2, j = Fl_add(j, i, p))
    1806       94491 :     k[j+1] = 1;
    1807        1764 :   a = n == 5 ? 0: k[1+Flx_lead(H)];
    1808        1764 :   d = Flx_difftable(H, p);
    1809      187577 :   for (i=0; i < p; i++)
    1810             :   {
    1811      185850 :     a += k[1+uel(d,n+1)];
    1812      185850 :     Flv_diffnext(d, p);
    1813             :   }
    1814        1727 :   set_avma(av);
    1815        1764 :   return a;
    1816             : }
    1817             : 
    1818             : static GEN
    1819        1889 : dirgenus2(GEN Q, GEN p, long n)
    1820             : {
    1821        1889 :   pari_sp av = avma;
    1822             :   GEN f;
    1823        1889 :   if (n > 2)
    1824         126 :     f = RgX_recip(hyperellcharpoly(gmul(Q,gmodulo(gen_1, p))));
    1825             :   else
    1826             :   {
    1827        1763 :     ulong pp = itou(p);
    1828        1763 :     GEN Qp = ZX_to_Flx(Q, pp);
    1829        1763 :     long t = Flx_genus2trace_naive(Qp, pp);
    1830        1764 :     f = deg1pol_shallow(stoi(t), gen_1, 0);
    1831             :   }
    1832        1890 :   return gerepileupto(av, RgXn_inv_i(f, n));
    1833             : }
    1834             : 
    1835             : GEN
    1836         756 : dirgenus2_worker(GEN P, ulong X, GEN Q)
    1837             : {
    1838         756 :   pari_sp av = avma;
    1839         756 :   long i, l = lg(P);
    1840         756 :   GEN V = cgetg(l, t_VEC);
    1841        2645 :   for(i = 1; i < l; i++)
    1842             :   {
    1843        1889 :     ulong p = uel(P,i);
    1844        1889 :     long d = ulogint(X, p) + 1; /* minimal d such that p^d > X */
    1845        1889 :     gel(V,i) = dirgenus2(Q, utoi(uel(P,i)), d);
    1846             :   }
    1847         756 :   return gerepilecopy(av, mkvec2(P,V));
    1848             : }
    1849             : 
    1850             : static GEN
    1851         168 : vecan_genus2(GEN an, long L)
    1852             : {
    1853         168 :   GEN Q = gel(an,1), bad = gel(an, 2);
    1854         168 :   GEN worker = snm_closure(is_entry("_dirgenus2_worker"), mkvec(Q));
    1855         168 :   return pardireuler(worker, gen_2, stoi(L), NULL, bad);
    1856             : }
    1857             : 
    1858             : static GEN
    1859          42 : genus2_redmodel(GEN P, GEN p)
    1860             : {
    1861          42 :   GEN M = FpX_factor(P, p);
    1862          42 :   GEN F = gel(M,1), E = gel(M,2);
    1863          42 :   long i, k, r = lg(F);
    1864          42 :   GEN U = scalarpol(leading_coeff(P), varn(P));
    1865          42 :   GEN G = cgetg(r, t_COL);
    1866         140 :   for (i=1, k=0; i<r; i++)
    1867             :   {
    1868          98 :     if (E[i]>1)
    1869          49 :       gel(G,++k) = gel(F,i);
    1870          98 :     if (odd(E[i]))
    1871          63 :       U = FpX_mul(U, gel(F,i), p);
    1872             :   }
    1873          42 :   setlg(G,++k);
    1874          42 :   return mkvec2(G,U);
    1875             : }
    1876             : 
    1877             : static GEN
    1878         280 : oneminusxd(long d)
    1879             : {
    1880         280 :   return gsub(gen_1, pol_xn(d, 0));
    1881             : }
    1882             : 
    1883             : static GEN
    1884          21 : ellfromeqncharpoly(GEN P, GEN Q, GEN p)
    1885             : {
    1886             :   long v;
    1887             :   GEN E, F, t, y;
    1888          21 :   v = fetch_var();
    1889          21 :   y = pol_x(v);
    1890          21 :   F = gsub(gadd(ZX_sqr(y), gmul(y, Q)), P);
    1891          21 :   E = ellinit(ellfromeqn(F), p, DEFAULTPREC);
    1892          21 :   delete_var();
    1893          21 :   t = ellap(E, p);
    1894          21 :   obj_free(E);
    1895          21 :   return mkpoln(3, gen_1, negi(t), p);
    1896             : }
    1897             : 
    1898             : static GEN
    1899          42 : genus2_eulerfact(GEN P, GEN p)
    1900             : {
    1901          42 :   GEN Pp = FpX_red(P, p);
    1902          42 :   GEN GU = genus2_redmodel(Pp, p);
    1903          42 :   long d = 6-degpol(Pp), v = d/2, w = odd(d);
    1904             :   GEN abe, tor;
    1905          42 :   GEN ki, kp = pol_1(0), kq = pol_1(0);
    1906          42 :   GEN F = gel(GU,1), Q = gel(GU,2);
    1907          42 :   long dQ = degpol(Q), lF = lg(F)-1;
    1908             : 
    1909           0 :   abe = dQ >= 5 ? RgX_recip(hyperellcharpoly(gmul(Q,gmodulo(gen_1,p))))
    1910          84 :       : dQ >= 3 ? RgX_recip(ellfromeqncharpoly(Q,gen_0,p))
    1911          42 :                 : pol_1(0);
    1912          35 :   ki = dQ != 0 ? oneminusxd(1)
    1913          49 :               : Fp_issquare(gel(Q,2),p) ? ZX_sqr(oneminusxd(1))
    1914           7 :                                         : oneminusxd(2);
    1915          42 :   if (lF)
    1916             :   {
    1917             :     long i;
    1918          84 :     for(i=1; i <= lF; i++)
    1919             :     {
    1920          49 :       GEN Fi = gel(F, i);
    1921          49 :       long d = degpol(Fi);
    1922          49 :       GEN e = FpX_rem(Q, Fi, p);
    1923          84 :       GEN kqf = lgpol(e)==0 ? oneminusxd(d):
    1924          70 :                 FpXQ_issquare(e, Fi, p) ? ZX_sqr(oneminusxd(d))
    1925          35 :                                         : oneminusxd(2*d);
    1926          49 :       kp = gmul(kp, oneminusxd(d));
    1927          49 :       kq = gmul(kq, kqf);
    1928             :     }
    1929             :   }
    1930          42 :   if (v)
    1931             :   {
    1932           7 :     GEN kqoo = w==1 ? oneminusxd(1):
    1933           0 :                Fp_issquare(leading_coeff(Q), p)? ZX_sqr(oneminusxd(1))
    1934           0 :                                               : oneminusxd(2);
    1935           7 :     kp = gmul(kp, oneminusxd(1));
    1936           7 :     kq = gmul(kq, kqoo);
    1937             :   }
    1938          42 :   tor = RgX_div(ZX_mul(oneminusxd(1), kq), ZX_mul(ki, kp));
    1939          42 :   return ginv( ZX_mul(abe, tor) );
    1940             : }
    1941             : 
    1942             : static GEN
    1943          28 : F2x_genus2_find_trans(GEN P, GEN Q, GEN F)
    1944             : {
    1945          28 :   pari_sp av = avma;
    1946          28 :   long i, d = F2x_degree(F), v = P[1];
    1947             :   GEN M, C, V;
    1948          28 :   M = cgetg(d+1, t_MAT);
    1949          84 :   for (i=1; i<=d; i++)
    1950             :   {
    1951          56 :     GEN Mi = F2x_rem(F2x_add(F2x_shift(Q,i-1), monomial_F2x(2*i-2,v)), F);
    1952          56 :     gel(M,i) = F2x_to_F2v(Mi, d);
    1953             :   }
    1954          28 :   C = F2x_to_F2v(F2x_rem(P, F), d);
    1955          28 :   V = F2m_F2c_invimage(M, C);
    1956          28 :   return gerepileuptoleaf(av, F2v_to_F2x(V, v));
    1957             : }
    1958             : 
    1959             : static GEN
    1960          35 : F2x_genus2_trans(GEN P, GEN Q, GEN H)
    1961             : {
    1962          35 :   return F2x_add(P,F2x_add(F2x_mul(H,Q), F2x_sqr(H)));
    1963             : }
    1964             : 
    1965             : static GEN
    1966          42 : F2x_genus_redoo(GEN P, GEN Q, long k)
    1967             : {
    1968          42 :   if (F2x_degree(P)==2*k)
    1969             :   {
    1970          14 :     long c = F2x_coeff(P,2*k-1), dQ = F2x_degree(Q);
    1971          14 :     if ((dQ==k-1 && c==1) || (dQ<k-1 && c==0))
    1972           7 :      return F2x_genus2_trans(P, Q, monomial_F2x(k, P[1]));
    1973             :   }
    1974          35 :   return P;
    1975             : }
    1976             : 
    1977             : static GEN
    1978          35 : F2x_pseudodisc(GEN P, GEN Q)
    1979             : {
    1980          35 :   GEN dP = F2x_deriv(P), dQ = F2x_deriv(Q);
    1981          35 :   return F2x_gcd(Q, F2x_add(F2x_mul(P, F2x_sqr(dQ)), F2x_sqr(dP)));
    1982             : }
    1983             : 
    1984             : static GEN
    1985          14 : F2x_genus_red(GEN P, GEN Q)
    1986             : {
    1987             :   long dP, dQ;
    1988             :   GEN F, FF;
    1989          14 :   P = F2x_genus_redoo(P, Q, 3);
    1990          14 :   P = F2x_genus_redoo(P, Q, 2);
    1991          14 :   P = F2x_genus_redoo(P, Q, 1);
    1992          14 :   dP = F2x_degree(P);
    1993          14 :   dQ = F2x_degree(Q);
    1994          14 :   FF = F = F2x_pseudodisc(P,Q);
    1995          35 :   while(F2x_degree(F)>0)
    1996             :   {
    1997          21 :     GEN M = gel(F2x_factor(F),1);
    1998          21 :     long i, l = lg(M);
    1999          49 :     for(i=1; i<l; i++)
    2000             :     {
    2001          28 :       GEN R = F2x_sqr(gel(M,i));
    2002          28 :       GEN H = F2x_genus2_find_trans(P, Q, R);
    2003          28 :       P = F2x_div(F2x_genus2_trans(P, Q, H), R);
    2004          28 :       Q = F2x_div(Q, gel(M,i));
    2005             :     }
    2006          21 :     F = F2x_pseudodisc(P, Q);
    2007             :   }
    2008          14 :   return mkvec4(P,Q,FF,mkvecsmall2(dP,dQ));
    2009             : }
    2010             : 
    2011             : /* Number of solutions of x^2+b*x+c */
    2012             : static long
    2013          21 : F2xqX_quad_nbroots(GEN b, GEN c, GEN T)
    2014             : {
    2015          21 :   if (lgpol(b) > 0)
    2016             :   {
    2017          14 :     GEN d = F2xq_div(c, F2xq_sqr(b, T), T);
    2018          14 :     return F2xq_trace(d, T)? 0: 2;
    2019             :   }
    2020             :   else
    2021           7 :     return 1;
    2022             : }
    2023             : 
    2024             : static GEN
    2025          28 : genus2_redmodel2(GEN P)
    2026             : {
    2027          28 :   GEN Q = pol_0(varn(P));
    2028          28 :   GEN P2 = ZX_to_F2x(P);
    2029          49 :   while (F2x_issquare(P2))
    2030             :   {
    2031          21 :     GEN H = F2x_to_ZX(F2x_sqrt(P2));
    2032          21 :     GEN P1 = ZX_sub(P, ZX_sqr(H));
    2033          21 :     GEN Q1 = ZX_add(Q, ZX_mulu(H, 2));
    2034          21 :     if ((signe(P1)==0 ||  ZX_lval(P1,2)>=2)
    2035          21 :      && (signe(Q1)==0 ||  ZX_lval(Q1,2)>=1))
    2036             :     {
    2037          21 :       P = ZX_shifti(P1, -2);
    2038          21 :       Q = ZX_shifti(Q1, -1);
    2039          21 :       P2= ZX_to_F2x(P);
    2040             :     } else break;
    2041             :   }
    2042          28 :   return mkvec2(P,Q);
    2043             : }
    2044             : 
    2045             : static GEN
    2046          14 : genus2_eulerfact2(GEN PQ)
    2047             : {
    2048          14 :   GEN V = F2x_genus_red(ZX_to_F2x(gel(PQ, 1)), ZX_to_F2x(gel(PQ, 2)));
    2049          14 :   GEN P = gel(V, 1), Q = gel(V, 2);
    2050          14 :   GEN F = gel(V, 3), v = gel(V, 4);
    2051             :   GEN abe, tor;
    2052          14 :   GEN ki, kp = pol_1(0), kq = pol_1(0);
    2053          14 :   long dP = F2x_degree(P), dQ = F2x_degree(Q), d = maxss(dP, 2*dQ);
    2054          14 :   if (!lgpol(F)) return pol_1(0);
    2055          21 :   ki = dQ!=0 || dP>0 ? oneminusxd(1):
    2056           7 :       dP==-1 ? ZX_sqr(oneminusxd(1)): oneminusxd(2);
    2057          28 :   abe = d>=5? RgX_recip(hyperellcharpoly(gmul(PQ,gmodulss(1,2)))):
    2058          14 :         d>=3? RgX_recip(ellfromeqncharpoly(F2x_to_ZX(P), F2x_to_ZX(Q), gen_2)):
    2059          14 :         pol_1(0);
    2060          14 :   if (lgpol(F))
    2061             :   {
    2062          14 :     GEN M = gel(F2x_factor(F), 1);
    2063          14 :     long i, lF = lg(M)-1;
    2064          35 :     for(i=1; i <= lF; i++)
    2065             :     {
    2066          21 :       GEN Fi = gel(M, i);
    2067          21 :       long d = F2x_degree(Fi);
    2068          21 :       long nb  = F2xqX_quad_nbroots(F2x_rem(Q, Fi), F2x_rem(P, Fi), Fi);
    2069          35 :       GEN kqf = nb==1 ? oneminusxd(d):
    2070           0 :                 nb==2 ? ZX_sqr(oneminusxd(d))
    2071          14 :                       : oneminusxd(2*d);
    2072          21 :       kp = gmul(kp, oneminusxd(d));
    2073          21 :       kq = gmul(kq, kqf);
    2074             :     }
    2075             :   }
    2076          14 :   if (maxss(v[1],2*v[2])<5)
    2077             :   {
    2078          14 :     GEN kqoo = v[1]>2*v[2] ? oneminusxd(1):
    2079           0 :                v[1]<2*v[2] ? ZX_sqr(oneminusxd(1))
    2080           7 :                            : oneminusxd(2);
    2081           7 :     kp = gmul(kp, oneminusxd(1));
    2082           7 :     kq = gmul(kq, kqoo);
    2083             :   }
    2084          14 :   tor = RgX_div(ZX_mul(oneminusxd(1),kq), ZX_mul(ki, kp));
    2085          14 :   return ginv( ZX_mul(abe, tor) );
    2086             : }
    2087             : 
    2088             : GEN
    2089          28 : lfungenus2(GEN G)
    2090             : {
    2091          28 :   pari_sp ltop = avma;
    2092             :   GEN Ldata;
    2093          28 :   GEN gr = genus2red(G, NULL);
    2094          28 :   GEN N  = gel(gr, 1), M = gel(gr, 2), Q = gel(gr, 3), L = gel(gr, 4);
    2095          28 :   GEN PQ = genus2_redmodel2(Q);
    2096             :   GEN e;
    2097          28 :   long i, lL = lg(L), ram2;
    2098          28 :   ram2 = absequaliu(gmael(M,1,1),2);
    2099          28 :   if (ram2 && equalis(gmael(M,2,1),-1))
    2100          14 :     pari_warn(warner,"unknown valuation of conductor at 2");
    2101          28 :   e = cgetg(lL+(ram2?0:1), t_VEC);
    2102          42 :   gel(e,1) = mkvec2(gen_2, ram2 ? genus2_eulerfact2(PQ)
    2103          14 :            : ginv( RgX_recip(hyperellcharpoly(gmul(PQ,gmodulss(1,2))))) );
    2104          70 :   for(i = ram2? 2: 1; i < lL; i++)
    2105             :   {
    2106          42 :     GEN Li = gel(L, i);
    2107          42 :     GEN p = gel(Li, 1);
    2108          42 :     gel(e, ram2 ? i: i+1) = mkvec2(p, genus2_eulerfact(Q,p));
    2109             :   }
    2110          28 :   Ldata = mkvecn(6, tag(mkvec2(Q,e), t_LFUN_GENUS2),
    2111             :       gen_0, mkvec4(gen_0, gen_0, gen_1, gen_1), gen_2, N, gen_0);
    2112          28 :   return gerepilecopy(ltop, Ldata);
    2113             : }
    2114             : 
    2115             : /*************************************************************/
    2116             : /*                        ETA QUOTIENTS                      */
    2117             : /* An eta quotient is a matrix with 2 columns [m, r_m] with  */
    2118             : /* m >= 1 representing f(\tau)=\prod_m\eta(m\tau)^{r_m}.     */
    2119             : /*************************************************************/
    2120             : 
    2121             : /* eta(x^v) + O(x^L) */
    2122             : GEN
    2123        1180 : eta_ZXn(long v, long L)
    2124             : {
    2125        1180 :   long n, k = 0, v2 = 2*v, bn = v, cn = 0;
    2126             :   GEN P;
    2127        1180 :   if (!L) return zeropol(0);
    2128        1180 :   P = cgetg(L+2,t_POL); P[1] = evalsigne(1);
    2129       64913 :   for(n = 0; n < L; n++) gel(P,n+2) = gen_0;
    2130        1180 :   for(n = 0;; n++, bn += v2, cn += v)
    2131        2762 :   { /* k = v * (3*n-1) * n / 2; bn = v * (2*n+1); cn = v * n */
    2132             :     long k2;
    2133        3942 :     gel(P, k+2) = odd(n)? gen_m1: gen_1;
    2134        3942 :     k2 = k+cn; if (k2 >= L) break;
    2135        3497 :     k = k2;
    2136             :     /* k = v * (3*n+1) * n / 2 */;
    2137        3497 :     gel(P, k+2) = odd(n)? gen_m1: gen_1;
    2138        3497 :     k2 = k+bn; if (k2 >= L) break;
    2139        2762 :     k = k2;
    2140             :   }
    2141        1180 :   setlg(P, k+3); return P;
    2142             : }
    2143             : GEN
    2144         308 : eta_product_ZXn(GEN eta, long L)
    2145             : {
    2146         308 :   pari_sp av = avma;
    2147         308 :   GEN P = NULL, D = gel(eta,1), R = gel(eta,2);
    2148         308 :   long i, l = lg(D);
    2149        1106 :   for (i = 1; i < l; ++i)
    2150             :   {
    2151         798 :     GEN Q = eta_ZXn(D[i], L);
    2152         798 :     long r = R[i];
    2153         798 :     if (r < 0) { Q = RgXn_inv_i(Q, L); r = -r; }
    2154         798 :     if (r != 1) Q = RgXn_powu_i(Q, r, L);
    2155         798 :     P = P? ZXn_mul(P, Q, L): Q;
    2156         798 :     if (gc_needed(av,1) && i > 1)
    2157             :     {
    2158           0 :       if (DEBUGMEM>1) pari_warn(warnmem,"eta_product_ZXn");
    2159           0 :       P = gerepilecopy(av, P);
    2160             :     }
    2161             :   }
    2162         308 :   return P;
    2163             : }
    2164             : static GEN
    2165         140 : vecan_eta(GEN an, long L)
    2166             : {
    2167         140 :   GEN v = RgX_to_RgC(eta_product_ZXn(an, L), L);
    2168         140 :   settyp(v, t_VEC); return v;
    2169             : }
    2170             : 
    2171             : /* return 1 if cuspidal, 0 if holomorphic, -1 otherwise */
    2172             : static int
    2173         210 : etacuspidal(GEN N, GEN k, GEN B, GEN R, GEN NB)
    2174             : {
    2175         210 :   long i, j, lD, l, cusp = 1;
    2176         210 :   pari_sp av = avma;
    2177             :   GEN D;
    2178         210 :   if (gsigne(k) < 0) return -1;
    2179         203 :   D = divisors(N); lD = lg(D); l = lg(B);
    2180        1246 :   for (i = 1; i < lD; i++)
    2181             :   {
    2182        1043 :     GEN t = gen_0, d = gel(D,i);
    2183             :     long s;
    2184        3283 :     for (j = 1; j < l; j++)
    2185        2240 :       t = addii(t, mulii(gel(NB,j), mulii(gel(R,j), sqri(gcdii(d, gel(B,j))))));
    2186        1043 :     s = signe(t);
    2187        1043 :     if (s < 0) return -1;
    2188        1043 :     if (!s) cusp = 0;
    2189             :   }
    2190         203 :   return gc_bool(av, cusp);
    2191             : }
    2192             : /* u | 24, level N = u*N0, N0 = lcm(B), NB[i] = N0/B[i] */
    2193             : static int
    2194          42 : etaselfdual(GEN B, GEN R, GEN NB, ulong u)
    2195             : {
    2196          42 :   pari_sp av = avma;
    2197          42 :   long i, l = lg(B);
    2198         140 :   for (i = 1; i < l; i++)
    2199             :   {
    2200          98 :     long j = ZV_search(B, muliu(gel(NB,i), u)); /* search for N / B[i] */
    2201          98 :     set_avma(av); if (!j || !equalii(gel(R,i),gel(R,j))) return 0;
    2202             :   }
    2203          42 :   return 1;
    2204             : }
    2205             : /* return Nebentypus of eta quotient, k2 = 2*k integral */
    2206             : static GEN
    2207         189 : etachar(GEN B, GEN R, GEN k2)
    2208             : {
    2209         189 :   long i, l = lg(B);
    2210         189 :   GEN P = gen_1;
    2211         504 :   for (i = 1; i < l; ++i) if (mpodd(gel(R,i))) P = mulii(P, gel(B,i));
    2212         189 :   switch(Mod4(k2))
    2213             :   {
    2214         119 :     case 0: break;
    2215          42 :     case 2:  P = negi(P); break;
    2216          28 :     default: P = shifti(P, 1); break;
    2217             :   }
    2218         189 :   return coredisc(P);
    2219             : }
    2220             : /* Return 0 if not on gamma_0(N). Sets conductor, modular weight, character,
    2221             :  * canonical matrix, v_q(eta), sd = 1 iff self-dual, cusp = 1 iff cuspidal
    2222             :  * [0 if holomorphic at all cusps, else -1] */
    2223             : long
    2224         238 : etaquotype(GEN *peta, GEN *pN, GEN *pk, GEN *CHI, long *pv, long *sd,
    2225             :            long *cusp)
    2226             : {
    2227         238 :   GEN B, R, S, T, N, NB, eta = *peta;
    2228             :   long l, i, u, S24;
    2229             : 
    2230         238 :   if (lg(eta) != 3) pari_err_TYPE("lfunetaquo", eta);
    2231         238 :   switch(typ(eta))
    2232             :   {
    2233          77 :     case t_VEC: eta = mkmat2(mkcol(gel(eta,1)), mkcol(gel(eta,2))); break;
    2234         161 :     case t_MAT: break;
    2235           0 :     default: pari_err_TYPE("lfunetaquo", eta);
    2236             :   }
    2237         238 :   if (!RgV_is_ZVpos(gel(eta,1)) || !RgV_is_ZV(gel(eta,2)))
    2238           0 :     pari_err_TYPE("lfunetaquo", eta);
    2239         238 :   *peta = eta = famat_reduce(eta);
    2240         238 :   B = gel(eta,1); l = lg(B); /* sorted in increasing order */
    2241         238 :   R = gel(eta,2);
    2242         238 :   N = ZV_lcm(B); NB = cgetg(l, t_VEC);
    2243         658 :   for (i = 1; i < l; i++) gel(NB,i) = diviiexact(N, gel(B,i));
    2244         238 :   S = gen_0; T = gen_0; u = 0;
    2245         658 :   for (i = 1; i < l; ++i)
    2246             :   {
    2247         420 :     GEN b = gel(B,i), r = gel(R,i);
    2248         420 :     S = addii(S, mulii(b, r));
    2249         420 :     T = addii(T, r);
    2250         420 :     u += umodiu(r,24) * umodiu(gel(NB,i), 24);
    2251             :   }
    2252         238 :   S = divis_rem(S, 24, &S24);
    2253         238 :   if (S24) return 0; /* non-integral valuation at oo */
    2254         231 :   u = 24 / ugcd(24, u % 24);
    2255         231 :   *pN = muliu(N, u); /* level */
    2256         231 :   *pk = gmul2n(T,-1); /* weight */
    2257         231 :   *pv = itos(S); /* valuation */
    2258         231 :   if (cusp) *cusp = etacuspidal(*pN, *pk, B, R, NB);
    2259         231 :   if (sd) *sd = etaselfdual(B, R, NB, u);
    2260         231 :   if (CHI) *CHI = etachar(B, R, T);
    2261         231 :   return 1;
    2262             : }
    2263             : 
    2264             : GEN
    2265          42 : lfunetaquo(GEN eta0)
    2266             : {
    2267          42 :   pari_sp ltop = avma;
    2268          42 :   GEN Ldata, N, BR, k, eta = eta0;
    2269             :   long v, sd, cusp;
    2270          42 :   if (!etaquotype(&eta, &N, &k, NULL, &v, &sd, &cusp))
    2271           0 :     pari_err_TYPE("lfunetaquo", eta0);
    2272          42 :   if (!cusp) pari_err_IMPL("noncuspidal eta quotient");
    2273          42 :   if (v != 1) pari_err_IMPL("valuation != 1");
    2274          42 :   if (!sd) pari_err_IMPL("non self-dual eta quotient");
    2275          42 :   if (typ(k) != t_INT) pari_err_TYPE("lfunetaquo [non-integral weight]", eta0);
    2276          42 :   BR = mkvec2(ZV_to_zv(gel(eta,1)), ZV_to_zv(gel(eta,2)));
    2277          42 :   Ldata = mkvecn(6, tag(BR,t_LFUN_ETA), gen_0, mkvec2(gen_0,gen_1), k,N, gen_1);
    2278          42 :   return gerepilecopy(ltop, Ldata);
    2279             : }
    2280             : 
    2281             : static GEN
    2282         455 : vecan_qf(GEN Q, long L)
    2283             : {
    2284         455 :   GEN v, w = qfrep0(Q, utoi(L), 1);
    2285             :   long i;
    2286         455 :   v = cgetg(L+1, t_VEC);
    2287       27524 :   for (i = 1; i <= L; i++) gel(v,i) = utoi(2 * w[i]);
    2288         455 :   return v;
    2289             : }
    2290             : 
    2291             : long
    2292         308 : qfiseven(GEN M)
    2293             : {
    2294         308 :   long i, l = lg(M);
    2295         735 :   for (i=1; i<l; i++)
    2296         644 :     if (mpodd(gcoeff(M,i,i))) return 0;
    2297          91 :   return 1;
    2298             : }
    2299             : 
    2300             : GEN
    2301          91 : lfunqf(GEN M, long prec)
    2302             : {
    2303          91 :   pari_sp ltop = avma;
    2304             :   long n;
    2305             :   GEN k, D, d, Mi, Ldata, poles, eno, dual;
    2306             : 
    2307          91 :   if (typ(M) != t_MAT) pari_err_TYPE("lfunqf", M);
    2308          91 :   if (!RgM_is_ZM(M))   pari_err_TYPE("lfunqf [not integral]", M);
    2309          91 :   n = lg(M)-1;
    2310          91 :   k = sstoQ(n,2);
    2311          91 :   M = Q_primpart(M);
    2312          91 :   Mi = ZM_inv(M, &d); /* d M^(-1) */
    2313          91 :   if (!qfiseven(M)) { M = gmul2n(M, 1); d = shifti(d,1); }
    2314          91 :   if (!qfiseven(Mi)){ Mi= gmul2n(Mi,1); d = shifti(d,1); }
    2315             :   /* det(Mi) = d^n/det(M), D^2 = det(Mi)/det(M) */
    2316          91 :   D = gdiv(gpow(d,k,prec), ZM_det(M));
    2317          91 :   if (!issquareall(D, &eno)) eno = gsqrt(D, prec);
    2318          91 :   dual = gequal1(D) ? gen_0: tag(Mi, t_LFUN_QF);
    2319          91 :   poles = mkcol2(mkvec2(k, simple_pole(gmul2n(eno,1))),
    2320             :                  mkvec2(gen_0, simple_pole(gen_m2)));
    2321          91 :   Ldata = mkvecn(7, tag(M, t_LFUN_QF), dual,
    2322             :        mkvec2(gen_0, gen_1), k, d, eno, poles);
    2323          91 :   return gerepilecopy(ltop, Ldata);
    2324             : }
    2325             : 
    2326             : /********************************************************************/
    2327             : /**  Artin L function, based on a GP script by Charlotte Euvrard   **/
    2328             : /********************************************************************/
    2329             : 
    2330             : static GEN
    2331         119 : artin_charfromgens(GEN G, GEN M)
    2332             : {
    2333         119 :   GEN R, V, ord = gal_get_orders(G), grp = gal_get_group(G);
    2334         119 :   long i, j, k, n = lg(ord)-1, m = lg(grp)-1;
    2335             : 
    2336         119 :   if (lg(M)-1 != n) pari_err_DIM("lfunartin");
    2337         119 :   R = cgetg(m+1, t_VEC);
    2338         119 :   gel(R, 1) = matid(lg(gel(M, 1))-1);
    2339         357 :   for (i = 1, k = 1; i <= n; ++i)
    2340             :   {
    2341         238 :     long c = k*(ord[i] - 1);
    2342         238 :     gel(R, ++k) = gel(M, i);
    2343        1043 :     for (j = 2; j <= c; ++j) gel(R, ++k) = gmul(gel(R,j), gel(M,i));
    2344             :   }
    2345         119 :   V = cgetg(m+1, t_VEC);
    2346        1281 :   for (i = 1; i <= m; i++) gel(V, gel(grp,i)[1]) = gtrace(gel(R,i));
    2347         119 :   return V;
    2348             : }
    2349             : 
    2350             : /* TODO move somewhere else? */
    2351             : GEN
    2352         140 : galois_get_conj(GEN G)
    2353             : {
    2354         140 :   GEN grp = gal_get_group(G);
    2355         140 :   long i, k, r = lg(grp)-1;
    2356         140 :   GEN b = zero_F2v(r);
    2357         392 :   for (k = 2; k <= r; ++k)
    2358             :   {
    2359         392 :     GEN g = gel(grp,k);
    2360         392 :     if (!F2v_coeff(b,g[1]) && g[g[1]]==1)
    2361             :     {
    2362         168 :       pari_sp av = avma;
    2363         168 :       GEN F = galoisfixedfield(G, g, 1, -1);
    2364         168 :       if (ZX_sturmpart(F, NULL) > 0) { set_avma(av); return g; }
    2365         140 :       for (i = 1; i<=r; i++)
    2366             :       {
    2367         112 :         GEN h = gel(grp, i);
    2368         112 :         long t = h[1];
    2369         112 :         while (h[t]!=1) t = h[t];
    2370         112 :         F2v_set(b, h[g[t]]);
    2371             :       }
    2372          28 :       set_avma(av);
    2373             :     }
    2374             :   }
    2375           0 :   pari_err_BUG("galois_get_conj");
    2376             :   return NULL;/*LCOV_EXCL_LINE*/
    2377             : }
    2378             : 
    2379        7700 : static GEN  cyclotoi(GEN v) { return simplify_shallow(lift_shallow(v)); }
    2380        2891 : static long cyclotos(GEN v) { return gtos(cyclotoi(v)); }
    2381        2821 : static long char_dim(GEN ch) { return cyclotos(gel(ch,1)); }
    2382             : 
    2383             : static GEN
    2384        1379 : artin_gamma(GEN N, GEN G, GEN ch)
    2385             : {
    2386        1379 :   long a, t, d = char_dim(ch);
    2387        1379 :   if (nf_get_r2(N) == 0) return vec01(d, 0);
    2388          70 :   a = galois_get_conj(G)[1];
    2389          70 :   t = cyclotos(gel(ch,a));
    2390          70 :   return vec01((d+t) / 2, (d-t) / 2);
    2391             : }
    2392             : 
    2393             : static long
    2394        3213 : artin_dim(GEN ind, GEN ch)
    2395             : {
    2396        3213 :   long n = lg(ch)-1;
    2397        3213 :   GEN elts = group_elts(ind, n);
    2398        3213 :   long i, d = lg(elts)-1;
    2399        3213 :   GEN s = gen_0;
    2400       18123 :   for(i=1; i<=d; i++)
    2401       14910 :     s = gadd(s, gel(ch, gel(elts,i)[1]));
    2402        3213 :   return gtos(gdivgs(cyclotoi(s), d));
    2403             : }
    2404             : 
    2405             : static GEN
    2406         623 : artin_ind(GEN elts, GEN ch, GEN p)
    2407             : {
    2408         623 :   long i, d = lg(elts)-1;
    2409         623 :   GEN s = gen_0;
    2410        2149 :   for(i=1; i<=d; i++)
    2411        1526 :     s = gadd(s, gel(ch, gmul(gel(elts,i),p)[1]));
    2412         623 :   return gdivgs(s, d);
    2413             : }
    2414             : 
    2415             : static GEN
    2416        2282 : artin_ram(GEN nf, GEN gal, GEN aut, GEN pr, GEN ramg, GEN ch, long d)
    2417             : {
    2418        2282 :   pari_sp av = avma;
    2419             :   long i, v, n;
    2420             :   GEN p, q, V, elts;
    2421        2282 :   if (d==0) return pol_1(0);
    2422         616 :   n = degpol(gal_get_pol(gal));
    2423         616 :   q = p = idealramfrobenius_aut(nf, gal, pr, ramg, aut);
    2424         616 :   elts = group_elts(gel(ramg,2), n);
    2425         616 :   v = fetch_var_higher();
    2426         616 :   V = cgetg(d+2, t_POL);
    2427         616 :   V[1] = evalsigne(1)|evalvarn(v);
    2428        1239 :   for(i=1; i<=d; i++)
    2429             :   {
    2430         623 :     gel(V,i+1) = artin_ind(elts, ch, q);
    2431         623 :     q = gmul(q, p);
    2432             :   }
    2433         616 :   delete_var();
    2434         616 :   V = RgXn_expint(RgX_neg(V),d+1);
    2435         616 :   setvarn(V,0); return gerepileupto(av, ginv(V));
    2436             : }
    2437             : 
    2438             : /* [Artin conductor, vec of [p, Lp]] */
    2439             : static GEN
    2440        1379 : artin_badprimes(GEN N, GEN G, GEN aut, GEN ch)
    2441             : {
    2442        1379 :   pari_sp av = avma;
    2443        1379 :   long i, d = char_dim(ch);
    2444        1379 :   GEN P = gel(absZ_factor(nf_get_disc(N)), 1);
    2445        1379 :   long lP = lg(P);
    2446        1379 :   GEN B = cgetg(lP, t_VEC), C = cgetg(lP, t_VEC);
    2447             : 
    2448        3661 :   for (i = 1; i < lP; ++i)
    2449             :   {
    2450        2282 :     GEN p = gel(P, i), pr = idealprimedec_galois(N, p);
    2451        2282 :     GEN J = idealramgroups_aut(N, G, pr, aut);
    2452        2282 :     GEN G0 = gel(J,2); /* inertia group */
    2453        2282 :     long lJ = lg(J);
    2454        2282 :     long sdec = artin_dim(G0, ch);
    2455        2282 :     long ndec = group_order(G0);
    2456        2282 :     long j, v = ndec * (d - sdec);
    2457        3213 :     for (j = 3; j < lJ; ++j)
    2458             :     {
    2459         931 :       GEN Jj = gel(J, j);
    2460         931 :       long s = artin_dim(Jj, ch);
    2461         931 :       v += group_order(Jj) * (d - s);
    2462             :     }
    2463        2282 :     gel(C, i) = powiu(p, v/ndec);
    2464        2282 :     gel(B, i) = mkvec2(p, artin_ram(N, G, aut, pr, J, ch, sdec));
    2465             :   }
    2466        1379 :   return gerepilecopy(av, mkvec2(ZV_prod(C), B));
    2467             : }
    2468             : 
    2469             : /* p does not divide nf.index */
    2470             : static GEN
    2471       52398 : idealfrobenius_easy(GEN nf, GEN gal, GEN aut, GEN T, GEN p)
    2472             : {
    2473       52398 :   long i, l = lg(aut), f = degpol(T);
    2474       52397 :   GEN D, Dzk, DzkT, DXp, grp = gal_get_group(gal);
    2475       52396 :   pari_sp av = avma;
    2476       52396 :   if (f==1) return gel(grp,1);
    2477       50583 :   Dzk = nf_get_zkprimpart(nf);
    2478       50581 :   D = modii(nf_get_zkden(nf), p);
    2479       50579 :   DzkT = RgV_to_RgM(FqV_red(Dzk, T, p), f);
    2480       50580 :   DXp = RgX_to_RgC(FpX_Frobenius(T, p), f);
    2481       50582 :   if (!equali1(D)) DXp = FpC_Fp_mul(DXp, D, p);
    2482      332076 :   for(i=1; i < l; i++)
    2483             :   {
    2484      332063 :     GEN g = gel(grp,i);
    2485      332063 :     if (perm_order(g)==f)
    2486             :     {
    2487      170227 :       GEN A = FpM_FpC_mul(DzkT, gel(aut,g[1]), p);
    2488      170202 :       if (ZV_equal(A, DXp)) {set_avma(av); return g; }
    2489             :     }
    2490             :   }
    2491             :   return NULL; /* LCOV_EXCL_LINE */
    2492             : }
    2493             : /* true nf; p divides nf.index, pr/p unramified */
    2494             : static GEN
    2495        1596 : idealfrobenius_hard(GEN nf, GEN gal, GEN aut, GEN pr)
    2496             : {
    2497        1596 :   long i, l = lg(aut), f = pr_get_f(pr);
    2498        1596 :   GEN modpr, p, T, X, Xp, pi, grp = gal_get_group(gal);
    2499        1596 :   pari_sp av = avma;
    2500        1596 :   if (f==1) return gel(grp,1);
    2501        1344 :   pi = pr_get_gen(pr);
    2502        1344 :   modpr = zkmodprinit(nf, pr);
    2503        1344 :   p = modpr_get_p(modpr);
    2504        1344 :   T = modpr_get_T(modpr);
    2505        1344 :   X = modpr_genFq(modpr);
    2506        1344 :   Xp = FpX_Frobenius(T, p);
    2507        9086 :   for (i = 1; i < l; i++)
    2508             :   {
    2509        9086 :     GEN g = gel(grp,i);
    2510        9086 :     if (perm_order(g)==f)
    2511             :     {
    2512        4270 :       GEN S = gel(aut,g[1]);
    2513        4270 :       GEN A = nf_to_Fq(nf, zk_galoisapplymod(nf,X,S,p), modpr);
    2514             :       /* sigma(X) = X^p (mod pr) and sigma(pi) in pr */
    2515        5726 :       if (ZX_equal(A, Xp) && (f == nf_get_degree(nf) ||
    2516        2800 :           ZC_prdvd(zk_galoisapplymod(nf,pi,S,p),pr))) { set_avma(av); return g; }
    2517             :     }
    2518             :   }
    2519             :   return NULL; /* LCOV_EXCL_LINE */
    2520             : }
    2521             : 
    2522             : static GEN
    2523       53993 : dirartin(GEN nf, GEN G, GEN V, GEN aut, GEN p, long n)
    2524             : {
    2525       53993 :   pari_sp av = avma;
    2526             :   GEN pr, frob;
    2527             :   /* pick one maximal ideal in the conjugacy class above p */
    2528       53993 :   GEN T = nf_get_pol(nf);
    2529       53992 :   if (!dvdii(nf_get_index(nf), p))
    2530             :   { /* simple case */
    2531       52394 :     GEN F = FpX_factor(T, p), P = gmael(F,1,1);
    2532       52399 :     frob = idealfrobenius_easy(nf, G, aut, P, p);
    2533             :   }
    2534             :   else
    2535             :   {
    2536        1596 :     pr = idealprimedec_galois(nf,p);
    2537        1596 :     frob = idealfrobenius_hard(nf, G, aut, pr);
    2538             :   }
    2539       53993 :   set_avma(av); return RgXn_inv(gel(V, frob[1]), n);
    2540             : }
    2541             : 
    2542             : GEN
    2543       15686 : dirartin_worker(GEN P, ulong X, GEN nf, GEN G, GEN V, GEN aut)
    2544             : {
    2545       15686 :   pari_sp av = avma;
    2546       15686 :   long i, l = lg(P);
    2547       15686 :   GEN W = cgetg(l, t_VEC);
    2548       69680 :   for(i = 1; i < l; i++)
    2549             :   {
    2550       53994 :     ulong p = uel(P,i);
    2551       53994 :     long d = ulogint(X, p) + 1; /* minimal d such that p^d > X */
    2552       53994 :     gel(W,i) = dirartin(nf, G, V, aut, utoi(uel(P,i)), d);
    2553             :   }
    2554       15686 :   return gerepilecopy(av, mkvec2(P,W));
    2555             : }
    2556             : 
    2557             : static GEN
    2558        2954 : vecan_artin(GEN an, long L, long prec)
    2559             : {
    2560        2954 :   GEN A, Sbad = gel(an,5);
    2561        2954 :   long n = itos(gel(an,6)), isreal = lg(an)<8 ? 0: !itos(gel(an,7));
    2562        2954 :   GEN worker = snm_closure(is_entry("_dirartin_worker"), vecslice(an,1,4));
    2563        2954 :   A = lift_shallow(pardireuler(worker, gen_2, stoi(L), NULL, Sbad));
    2564        2954 :   A = RgXV_RgV_eval(A, grootsof1(n, prec));
    2565        2954 :   if (isreal) A = real_i(A);
    2566        2954 :   return A;
    2567             : }
    2568             : 
    2569             : static GEN
    2570        2856 : char_expand(GEN conj, GEN ch)
    2571             : {
    2572        2856 :   long i, l = lg(conj);
    2573        2856 :   GEN V = cgetg(l, t_COL);
    2574       31409 :   for (i=1; i<l; i++) gel(V,i) = gel(ch, conj[i]);
    2575        2856 :   return V;
    2576             : }
    2577             : 
    2578             : static GEN
    2579        1596 : handle_zeta(long n, GEN ch, long *m)
    2580             : {
    2581             :   GEN c;
    2582        1596 :   long t, i, l = lg(ch);
    2583        1596 :   GEN dim = cyclotoi(vecsum(ch));
    2584        1596 :   if (typ(dim) != t_INT)
    2585           0 :     pari_err_DOMAIN("lfunartin","chi","is not a", strtoGENstr("character"), ch);
    2586        1596 :   t = itos(dim);
    2587        1596 :   if (t < 0 || t % n)
    2588           0 :     pari_err_DOMAIN("lfunartin","chi","is not a", strtoGENstr("character"), ch);
    2589        1596 :   if (t == 0) { *m = 0; return ch; }
    2590         224 :   *m = t / n;
    2591         224 :   c = cgetg(l, t_COL);
    2592        2065 :   for (i=1; i<l; i++)
    2593        1841 :     gel(c,i) = gsubgs(gel(ch,i), *m);
    2594         224 :   return c;
    2595             : }
    2596             : 
    2597             : static int
    2598        6496 : cyclo_is_real(GEN v, GEN ix)
    2599             : {
    2600        6496 :   pari_sp av = avma;
    2601        6496 :   GEN w = poleval(lift_shallow(v), ix);
    2602        6496 :   return gc_bool(av, gequal(w, v));
    2603             : }
    2604             : 
    2605             : static int
    2606        1379 : char_is_real(GEN ch, GEN mod)
    2607             : {
    2608        1379 :   long i, l = lg(ch);
    2609        1379 :   GEN ix = QXQ_inv(pol_x(varn(mod)), mod);
    2610        7014 :   for (i=1; i<l; i++)
    2611        6496 :     if (!cyclo_is_real(gel(ch,i), ix)) return 0;
    2612         518 :   return 1;
    2613             : }
    2614             : 
    2615             : GEN
    2616        1610 : lfunartin(GEN nf, GEN gal, GEN ch, long o, long bitprec)
    2617             : {
    2618        1610 :   pari_sp av = avma;
    2619        1610 :   GEN bc, V, aut, mod, Ldata = NULL, chx, cc, conj, repr;
    2620             :   long tmult, var;
    2621        1610 :   nf = checknf(nf);
    2622        1610 :   checkgal(gal);
    2623        1610 :   var = gvar(ch);
    2624        1610 :   if (var == 0) pari_err_PRIORITY("lfunartin",ch,"=",0);
    2625        1610 :   if (var < 0) var = 1;
    2626        1610 :   if (!is_vec_t(typ(ch))) pari_err_TYPE("lfunartin", ch);
    2627        1610 :   cc = group_to_cc(gal);
    2628        1610 :   conj = gel(cc,2);
    2629        1610 :   repr = gel(cc,3);
    2630        1610 :   mod = mkpolmod(gen_1, polcyclo(o, var));
    2631        1610 :   if (lg(ch)>1 && typ(gel(ch,1))==t_MAT)
    2632         119 :     chx = artin_charfromgens(gal, gmul(ch,mod));
    2633             :   else
    2634             :   {
    2635        1491 :     if (lg(repr) != lg(ch)) pari_err_DIM("lfunartin");
    2636        1477 :     chx = char_expand(conj, gmul(ch,mod));
    2637             :   }
    2638        1596 :   chx = handle_zeta(nf_get_degree(nf), chx, &tmult);
    2639        1596 :   ch = shallowextract(chx, repr);
    2640        1596 :   if (!gequal0(chx))
    2641             :   {
    2642        1379 :     GEN real = char_is_real(chx, gel(mod,1))? gen_0: gen_1;
    2643        1379 :     aut = nfgaloispermtobasis(nf, gal);
    2644        1379 :     V = gmul(char_expand(conj, galoischarpoly(gal, ch, o)), mod);
    2645        1379 :     bc = artin_badprimes(nf, gal, aut, chx);
    2646        2758 :     Ldata = mkvecn(6,
    2647        1379 :       tag(mkcoln(7, nf, gal, V, aut, gel(bc, 2), stoi(o), real), t_LFUN_ARTIN),
    2648        1379 :       real, artin_gamma(nf, gal, chx), gen_1, gel(bc,1), gen_0);
    2649             :   }
    2650        1596 :   if (tmult==0 && Ldata==NULL) pari_err_TYPE("lfunartin",ch);
    2651        1596 :   if (tmult)
    2652             :   {
    2653             :     long i;
    2654         224 :     if (Ldata==NULL) { Ldata = lfunzeta(); tmult--; }
    2655         231 :     for(i=1; i<=tmult; i++)
    2656           7 :       Ldata = lfunmul(Ldata, gen_1, bitprec);
    2657             :   }
    2658        1596 :   return gerepilecopy(av, Ldata);
    2659             : }
    2660             : 
    2661             : static GEN
    2662          21 : lfunzetakinit_artin(GEN nf, GEN gal, GEN dom, long der, long bitprec)
    2663             : {
    2664          21 :   pari_sp ltop = avma;
    2665          21 :   GEN To = galoischartable(gal), T = gel(To, 1);
    2666          21 :   long o = itos(gel(To, 2));
    2667             :   GEN F, E, M, domain;
    2668          21 :   long i, l = lg(T);
    2669          21 :   F = cgetg(l, t_VEC);
    2670          21 :   E = cgetg(l, t_VECSMALL);
    2671          84 :   for (i = 1; i < l; ++i)
    2672             :   {
    2673          63 :     GEN L = lfunartin(nf, gal, gel(T,i), o, bitprec);
    2674          63 :     gel(F, i) = lfuninit(L, dom, der, bitprec);
    2675          63 :     E[i] = char_dim(gel(T,i));
    2676             :   }
    2677          21 :   domain = mkvec2(dom, mkvecsmall2(der, bitprec));
    2678          21 :   M = mkvec3(F, E, zero_zv(l-1));
    2679          21 :   return gerepilecopy(ltop, lfuninit_make(t_LDESC_PRODUCT, lfunzetak_i(nf),
    2680             :                                           M, domain));
    2681             : }
    2682             : 
    2683             : /********************************************************************/
    2684             : /**                    High-level Constructors                     **/
    2685             : /********************************************************************/
    2686             : enum { t_LFUNMISC_POL, t_LFUNMISC_CHIQUAD, t_LFUNMISC_CHICONREY,
    2687             :        t_LFUNMISC_CHIGEN, t_LFUNMISC_ELLINIT, t_LFUNMISC_ETAQUO };
    2688             : static long
    2689        8036 : lfundatatype(GEN data)
    2690             : {
    2691        8036 :   switch(typ(data))
    2692             :   {
    2693        3360 :     case t_INT: return t_LFUNMISC_CHIQUAD;
    2694         140 :     case t_INTMOD: return t_LFUNMISC_CHICONREY;
    2695         497 :     case t_POL: return t_LFUNMISC_POL;
    2696        4039 :     case t_VEC:
    2697        4039 :       if (checknf_i(data)) return t_LFUNMISC_POL;
    2698        3920 :       switch(lg(data))
    2699             :       {
    2700        1736 :         case 17: return t_LFUNMISC_ELLINIT;
    2701        1883 :         case 3:  return t_LFUNMISC_CHIGEN;
    2702             :       }
    2703         301 :       break;
    2704             :   }
    2705         301 :   return -1;
    2706             : }
    2707             : static GEN
    2708       64057 : lfunmisc_to_ldata_i(GEN ldata, long shallow)
    2709             : {
    2710             :   long lx;
    2711       64057 :   if (is_linit(ldata)) ldata = linit_get_ldata(ldata);
    2712       64057 :   lx = lg(ldata);
    2713       64057 :   if (typ(ldata)==t_VEC && (lx == 7 || lx == 8) && is_tagged(ldata))
    2714             :   {
    2715       56021 :     if (!shallow) ldata = gcopy(ldata);
    2716       56021 :     checkldata(ldata); return ldata;
    2717             :   }
    2718        8036 :   switch (lfundatatype(ldata))
    2719             :   {
    2720         616 :     case t_LFUNMISC_POL: return lfunzetak(ldata);
    2721        3360 :     case t_LFUNMISC_CHIQUAD: return lfunchiquad(ldata);
    2722         140 :     case t_LFUNMISC_CHICONREY:
    2723             :     {
    2724         140 :       GEN G = znstar0(gel(ldata,1), 1);
    2725         140 :       return lfunchiZ(G, gel(ldata,2));
    2726             :     }
    2727        1883 :     case t_LFUNMISC_CHIGEN:
    2728             :     {
    2729        1883 :       GEN G = gel(ldata,1), chi = gel(ldata,2);
    2730        1883 :       switch(nftyp(G))
    2731             :       {
    2732        1582 :         case typ_BIDZ: return lfunchiZ(G, chi);
    2733         133 :         case typ_BNR: return lfunchigen(G, chi);
    2734             :       }
    2735             :     }
    2736         168 :     break;
    2737             : 
    2738        1736 :     case t_LFUNMISC_ELLINIT: return lfunell(ldata);
    2739             :   }
    2740         469 :   if (shallow != 2) pari_err_TYPE("lfunmisc_to_ldata",ldata);
    2741         462 :   return NULL;
    2742             : }
    2743             : 
    2744             : GEN
    2745        1057 : lfunmisc_to_ldata(GEN ldata)
    2746        1057 : { return lfunmisc_to_ldata_i(ldata, 0); }
    2747             : 
    2748             : GEN
    2749       49840 : lfunmisc_to_ldata_shallow(GEN ldata)
    2750       49840 : { return lfunmisc_to_ldata_i(ldata, 1); }
    2751             : 
    2752             : GEN
    2753       13160 : lfunmisc_to_ldata_shallow_i(GEN ldata)
    2754       13160 : { return lfunmisc_to_ldata_i(ldata, 2); }
    2755             : 
    2756             : /********************************************************************/
    2757             : /**                    High-level an expansion                     **/
    2758             : /********************************************************************/
    2759             : /* van is the output of ldata_get_an: return a_1,...a_L at precision prec */
    2760             : GEN
    2761       17801 : ldata_vecan(GEN van, long L, long prec)
    2762             : {
    2763       17801 :   GEN an = gel(van, 2);
    2764       17801 :   long t = mael(van,1,1);
    2765             :   pari_timer ti;
    2766       17801 :   if (DEBUGLEVEL >= 1)
    2767           0 :     err_printf("Lfun: computing %ld coeffs, prec %ld, type %ld\n", L, prec, t);
    2768       17801 :   if (DEBUGLEVEL >= 2) timer_start(&ti);
    2769       17801 :   switch (t)
    2770             :   {
    2771             :     long n;
    2772        1890 :     case t_LFUN_GENERIC:
    2773        1890 :       an = vecan_closure(an, L, prec);
    2774        1820 :       n = lg(an)-1;
    2775        1820 :       if (n < L)
    2776             :       {
    2777          14 :         pari_warn(warner, "#an = %ld < %ld, results may be imprecise", n, L);
    2778          14 :         an = shallowconcat(an, zerovec(L-n));
    2779             :       }
    2780        1820 :       break;
    2781           0 :     case t_LFUN_CLOSURE0:
    2782             :       pari_err_BUG("ldata_vecan: please call ldata_newprec");/*LCOV_EXCL_LINE*/
    2783        1904 :     case t_LFUN_ZETA: an = const_vecsmall(L, 1); break;
    2784        2156 :     case t_LFUN_NF:  an = dirzetak(an, stoi(L)); break;
    2785        1967 :     case t_LFUN_ELL:
    2786        1967 :       an = (ell_get_type(an) == t_ELL_Q) ? ellanQ_zv(an, L): ellan(an, L);
    2787        1967 :       break;
    2788        1288 :     case t_LFUN_KRONECKER: an = vecan_Kronecker(an, L); break;
    2789         987 :     case t_LFUN_CHIZ: an = vecan_chiZ(an, L, prec); break;
    2790         910 :     case t_LFUN_CHIGEN: an = vecan_chigen(an, L, prec); break;
    2791        2954 :     case t_LFUN_ARTIN: an = vecan_artin(an, L, prec); break;
    2792         140 :     case t_LFUN_ETA: an = vecan_eta(an, L); break;
    2793         455 :     case t_LFUN_QF: an = vecan_qf(an, L); break;
    2794         637 :     case t_LFUN_DIV: an = vecan_div(an, L, prec); break;
    2795         315 :     case t_LFUN_MUL: an = vecan_mul(an, L, prec); break;
    2796          98 :     case t_LFUN_CONJ: an = vecan_conj(an, L, prec); break;
    2797         343 :     case t_LFUN_SYMPOW_ELL: an = vecan_ellsympow(an, L); break;
    2798         168 :     case t_LFUN_GENUS2: an = vecan_genus2(an, L); break;
    2799         406 :     case t_LFUN_MFCLOS:
    2800             :     {
    2801         406 :       GEN F = gel(an,1), E = gel(an,2), c = gel(an,3);
    2802         406 :       an = mfcoefs(F,L,1) + 1; /* skip a_0 */
    2803         406 :       an[0] = evaltyp(t_VEC)|evallg(L+1);
    2804         406 :       an = mfvecembed(E, an);
    2805         406 :       if (!isint1(c)) an = RgV_Rg_mul(an,c);
    2806         406 :       break;
    2807             :     }
    2808         581 :     case t_LFUN_TWIST: an = vecan_twist(an, L, prec); break;
    2809         602 :     case t_LFUN_SHIFT: an = vecan_shift(an, L, prec); break;
    2810           0 :     default: pari_err_TYPE("ldata_vecan", van);
    2811             :   }
    2812       17731 :   if (DEBUGLEVEL >= 2) timer_printf(&ti, "ldata_vecan");
    2813       17731 :   return an;
    2814             : }
    2815             : 
    2816             : /* shallow function */
    2817             : GEN
    2818       17535 : ldata_newprec(GEN ldata, long prec)
    2819             : {
    2820       17535 :   GEN van = ldata_get_an(ldata);
    2821       17535 :   GEN an = gel(van, 2);
    2822       17535 :   long t = mael(van,1,1);
    2823       17535 :   switch (t)
    2824             :   {
    2825         140 :     case t_LFUN_CLOSURE0:
    2826         140 :       return lfuncreate(closure_callgen0prec(an, prec));
    2827         329 :     case t_LFUN_QF:
    2828             :     {
    2829         329 :       GEN eno = ldata_get_rootno(ldata);
    2830         329 :       if (typ(eno)==t_REAL && realprec(eno) < prec)
    2831          56 :         return lfunqf(an, prec);
    2832         273 :       break;
    2833             :     }
    2834             :   }
    2835       17339 :   return ldata;
    2836             : }

Generated by: LCOV version 1.13