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.14.0 lcov report (development 27775-aca467eab2) Lines: 1746 1893 92.2 %
Date: 2022-07-03 07:33:15 Functions: 168 178 94.4 %
Legend: Lines: hit not hit

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

Generated by: LCOV version 1.13