Code coverage tests

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

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

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

LCOV - code coverage report
Current view: top level - basemath - lfunutils.c (source / functions) Hit Total Coverage
Test: PARI/GP v2.10.0 lcov report (development 21061-8feaff2) Lines: 1155 1245 92.8 %
Date: 2017-09-24 06:24:57 Functions: 104 108 96.3 %
Legend: Lines: hit not hit

          Line data    Source code
       1             : /* Copyright (C) 2015  The PARI group.
       2             : 
       3             : This file is part of the PARI/GP package.
       4             : 
       5             : PARI/GP is free software; you can redistribute it and/or modify it under the
       6             : terms of the GNU General Public License as published by the Free Software
       7             : Foundation. It is distributed in the hope that it will be useful, but WITHOUT
       8             : ANY WARRANTY WHATSOEVER.
       9             : 
      10             : Check the License for details. You should have received a copy of it, along
      11             : with the package; see the file 'COPYING'. If not, write to the Free Software
      12             : Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. */
      13             : 
      14             : /********************************************************************/
      15             : /**                                                                **/
      16             : /**                 L-functions: Applications                      **/
      17             : /**                                                                **/
      18             : /********************************************************************/
      19             : 
      20             : #include "pari.h"
      21             : #include "paripriv.h"
      22             : 
      23             : static GEN
      24        6993 : tag(GEN x, long t) { return mkvec2(mkvecsmall(t), x); }
      25             : 
      26             : /* v a t_VEC of length > 1 */
      27             : static int
      28       19479 : is_tagged(GEN v)
      29             : {
      30       19479 :   GEN T = gel(v,1);
      31       19479 :   return (typ(T)==t_VEC && lg(T)==3 && typ(gel(T,1))==t_VECSMALL);
      32             : }
      33             : static void
      34       19479 : checkldata(GEN ldata)
      35             : {
      36             :   GEN vga, w, N;
      37             : #if 0 /* assumed already checked and true */
      38             :   long l = lg(ldata);
      39             :   if (typ(ldata)!=t_VEC || l < 7 || l > 8 || !is_tagged(ldata))
      40             :     pari_err_TYPE("checkldata", ldata);
      41             : #endif
      42       19479 :   vga = ldata_get_gammavec(ldata);
      43       19479 :   if (typ(vga) != t_VEC) pari_err_TYPE("checkldata [gammavec]",vga);
      44       19479 :   w = gel(ldata, 4); /* FIXME */
      45       19479 :   switch(typ(w))
      46             :   {
      47       19465 :     case t_INT: break;
      48          14 :     case t_VEC: if (lg(w) == 3 && typ(gel(w,1)) == t_INT) break;
      49           0 :     default: pari_err_TYPE("checkldata [weight]",w);
      50             :   }
      51       19479 :   N = ldata_get_conductor(ldata);
      52       19479 :   if (typ(N) != t_INT) pari_err_TYPE("checkldata [conductor]",N);
      53       19479 : }
      54             : 
      55             : /* data may be either an object (polynomial, elliptic curve, etc...)
      56             :  * or a description vector [an,sd,Vga,k,conductor,rootno,{poles}]. */
      57             : GEN
      58         539 : lfuncreate(GEN data)
      59             : {
      60         539 :   long lx = lg(data);
      61         539 :   if (typ(data)==t_VEC && (lx == 7 || lx == 8))
      62             :   {
      63             :     GEN ldata;
      64         322 :     if (is_tagged(data)) ldata = gcopy(data);
      65             :     else
      66             :     { /* tag first component as t_LFUN_GENERIC */
      67         252 :       ldata = gcopy(data);
      68         252 :       gel(ldata, 1) = tag(gel(ldata,1), t_LFUN_GENERIC);
      69         252 :       if (typ(gel(ldata, 2))!=t_INT)
      70          21 :         gel(ldata, 2) = tag(gel(ldata,2), t_LFUN_GENERIC);
      71             :     }
      72         322 :     checkldata(ldata); return ldata;
      73             :   }
      74         217 :   return lfunmisc_to_ldata(data);
      75             : }
      76             : 
      77             : /********************************************************************/
      78             : /**                     Simple constructors                        **/
      79             : /********************************************************************/
      80             : 
      81             : static GEN
      82           0 : vecan_conj(GEN an, long n, long prec)
      83             : {
      84           0 :   GEN p1 = ldata_vecan(gel(an,1), n, prec);
      85           0 :   return typ(p1) == t_VEC? gconj(p1): p1;
      86             : }
      87             : 
      88             : static GEN
      89         161 : vecan_mul(GEN an, long n, long prec)
      90             : {
      91         161 :   GEN p1 = ldata_vecan(gel(an,1), n, prec);
      92         161 :   GEN p2 = ldata_vecan(gel(an,2), n, prec);
      93         161 :   if (typ(p1) == t_VECSMALL) p1 = vecsmall_to_vec(p1);
      94         161 :   if (typ(p2) == t_VECSMALL) p2 = vecsmall_to_vec(p2);
      95         161 :   return dirmul(p1, p2);
      96             : }
      97             : 
      98             : static GEN
      99          35 : lfunconvol(GEN a1, GEN a2)
     100          35 : { return tag(mkvec2(a1, a2), t_LFUN_MUL); }
     101             : 
     102             : static GEN
     103         609 : vecan_div(GEN an, long n, long prec)
     104             : {
     105         609 :   GEN p1 = ldata_vecan(gel(an,1), n, prec);
     106         609 :   GEN p2 = ldata_vecan(gel(an,2), n, prec);
     107         609 :   if (typ(p1) == t_VECSMALL) p1 = vecsmall_to_vec(p1);
     108         609 :   if (typ(p2) == t_VECSMALL) p2 = vecsmall_to_vec(p2);
     109         609 :   return dirdiv(p1, p2);
     110             : }
     111             : 
     112             : static GEN
     113          49 : lfunconvolinv(GEN a1, GEN a2)
     114          49 : { return tag(mkvec2(a1,a2), t_LFUN_DIV); }
     115             : 
     116             : static GEN
     117           0 : lfunconj(GEN a1)
     118           0 : { return tag(mkvec(a1), t_LFUN_CONJ); }
     119             : 
     120             : static GEN
     121          84 : lfuncombdual(GEN fun(GEN, GEN), GEN ldata1, GEN ldata2)
     122             : {
     123          84 :   GEN a1 = ldata_get_an(ldata1), a2 = ldata_get_an(ldata2);
     124          84 :   GEN b1 = ldata_get_dual(ldata1), b2 = ldata_get_dual(ldata2);
     125          84 :   if (typ(b1)==t_INT && typ(b2)==t_INT)
     126          84 :     return utoi(signe(b1) && signe(b2));
     127             :   else
     128             :   {
     129           0 :     if (typ(b1)==t_INT) b1 = signe(b1) ? lfunconj(a1): a1;
     130           0 :     if (typ(b2)==t_INT) b2 = signe(b2) ? lfunconj(a2): a2;
     131           0 :     return fun(b1, b2);
     132             :   }
     133             : }
     134             : 
     135             : static GEN
     136          35 : lfunmulpoles(GEN ldata1, GEN ldata2, long bitprec)
     137             : {
     138          35 :   long k = ldata_get_k(ldata1), l, j;
     139          35 :   GEN r1 = ldata_get_residue(ldata1);
     140          35 :   GEN r2 = ldata_get_residue(ldata2), r;
     141             : 
     142          35 :   if (r1 && typ(r1) != t_VEC) r1 = mkvec(mkvec2(stoi(k), r1));
     143          35 :   if (r2 && typ(r2) != t_VEC) r2 = mkvec(mkvec2(stoi(k), r2));
     144          35 :   if (!r1)
     145             :   {
     146           7 :     if (!r2) return NULL;
     147           0 :     r1 = lfunrtopoles(r2);
     148             :   }
     149             :   else
     150             :   {
     151          28 :     r1 = lfunrtopoles(r1);
     152          28 :     if (r2) r1 = setunion(r1, lfunrtopoles(r2));
     153             :   }
     154          28 :   l = lg(r1); r = cgetg(l, t_VEC);
     155          56 :   for (j = 1; j < l; j++)
     156             :   {
     157          28 :     GEN be = gel(r1,j);
     158          28 :     GEN z1 = lfun(ldata1,be,bitprec), z2 = lfun(ldata2,be,bitprec);
     159          28 :     if (typ(z1) == t_SER && typ(z2) == t_SER)
     160             :     { /* pole of both, recompute to needed seriesprecision */
     161          28 :       long e = valp(z1) + valp(z2);
     162          28 :       GEN b = RgX_to_ser(deg1pol_shallow(gen_1, be, 0), 3-e);
     163          28 :       z1 = lfun(ldata1,b,bitprec);
     164          28 :       z2 = lfun(ldata2,b,bitprec);
     165             :     }
     166          28 :     gel(r,j) = mkvec2(be, gmul(z1, z2));
     167             :   }
     168          28 :   return r;
     169             : }
     170             : 
     171             : GEN
     172          35 : lfunmul(GEN ldata1, GEN ldata2, long bitprec)
     173             : {
     174          35 :   pari_sp ltop = avma;
     175             :   GEN r, N, Vga, eno, a1a2, b1b2, LD;
     176             :   long k;
     177          35 :   ldata1 = lfunmisc_to_ldata_shallow(ldata1);
     178          35 :   ldata2 = lfunmisc_to_ldata_shallow(ldata2);
     179          35 :   k = ldata_get_k(ldata1);
     180          35 :   if (ldata_get_k(ldata2) != k)
     181           0 :     pari_err_OP("lfunmul [weight]",ldata1, ldata2);
     182          35 :   r = lfunmulpoles(ldata1, ldata2, bitprec);
     183          35 :   N = gmul(ldata_get_conductor(ldata1), ldata_get_conductor(ldata2));
     184          35 :   Vga = shallowconcat(ldata_get_gammavec(ldata1), ldata_get_gammavec(ldata2));
     185          35 :   Vga = sort(Vga);
     186          35 :   eno = gmul(ldata_get_rootno(ldata1), ldata_get_rootno(ldata2));
     187          35 :   a1a2 = lfunconvol(ldata_get_an(ldata1), ldata_get_an(ldata2));
     188          35 :   b1b2 = lfuncombdual(lfunconvol, ldata1, ldata2);
     189          35 :   LD = mkvecn(7, a1a2, b1b2, Vga, stoi(k), N, eno, r);
     190          35 :   if (!r) setlg(LD,7);
     191          35 :   return gerepilecopy(ltop, LD);
     192             : }
     193             : 
     194             : static GEN
     195          49 : lfundivpoles(GEN ldata1, GEN ldata2, long bitprec)
     196             : {
     197          49 :   long k = ldata_get_k(ldata1), i, j, l;
     198          49 :   GEN r1 = ldata_get_residue(ldata1);
     199          49 :   GEN r2 = ldata_get_residue(ldata2), r;
     200             : 
     201          49 :   if (r1 && typ(r1) != t_VEC) r1 = mkvec(mkvec2(stoi(k), r1));
     202          49 :   if (r2 && typ(r2) != t_VEC) r2 = mkvec(mkvec2(stoi(k), r2));
     203          49 :   if (!r1) return NULL;
     204          49 :   r1 = lfunrtopoles(r1);
     205          49 :   l = lg(r1); r = cgetg(l, t_VEC);
     206          98 :   for (i = j = 1; j < l; j++)
     207             :   {
     208          49 :     GEN be = gel(r1,j);
     209          49 :     GEN z = gdiv(lfun(ldata1,be,bitprec), lfun(ldata2,be,bitprec));
     210          49 :     if (valp(z) < 0) gel(r,i++) = mkvec2(be, z);
     211             :   }
     212          49 :   if (i == 1) return NULL;
     213          14 :   setlg(r, i); return r;
     214             : }
     215             : 
     216             : GEN
     217          49 : lfundiv(GEN ldata1, GEN ldata2, long bitprec)
     218             : {
     219          49 :   pari_sp ltop = avma;
     220             :   GEN r, N, v, v1, v2, eno, a1a2, b1b2, LD, eno2;
     221             :   long k, j, j1, j2, l1, l2;
     222          49 :   ldata1 = lfunmisc_to_ldata_shallow(ldata1);
     223          49 :   ldata2 = lfunmisc_to_ldata_shallow(ldata2);
     224          49 :   k = ldata_get_k(ldata1);
     225          49 :   if (ldata_get_k(ldata2) != k)
     226           0 :     pari_err_OP("lfundiv [weight]",ldata1, ldata2);
     227          49 :   r = lfundivpoles(ldata1, ldata2, bitprec);
     228          49 :   N = gdiv(ldata_get_conductor(ldata1), ldata_get_conductor(ldata2));
     229          49 :   if (typ(N) != t_INT) pari_err_OP("lfundiv [conductor]",ldata1, ldata2);
     230          49 :   a1a2 = lfunconvolinv(ldata_get_an(ldata1), ldata_get_an(ldata2));
     231          49 :   b1b2 = lfuncombdual(lfunconvolinv, ldata1, ldata2);
     232          49 :   eno2 = ldata_get_rootno(ldata2);
     233          49 :   eno = isintzero(eno2)? gen_0: gdiv(ldata_get_rootno(ldata1), eno2);
     234          49 :   v1 = shallowcopy(ldata_get_gammavec(ldata1));
     235          49 :   v2 = ldata_get_gammavec(ldata2);
     236          49 :   l1 = lg(v1); l2 = lg(v2);
     237         105 :   for (j2 = 1; j2 < l2; j2++)
     238             :   {
     239          84 :     for (j1 = 1; j1 < l1; j1++)
     240          84 :       if (gel(v1,j1) && gequal(gel(v1,j1), gel(v2,j2)))
     241             :       {
     242          56 :         gel(v1,j1) = NULL; break;
     243             :       }
     244          56 :     if (j1 == l1) pari_err_OP("lfundiv [Vga]",ldata1, ldata2);
     245             :   }
     246          49 :   v = cgetg(l1-l2+1, t_VEC);
     247         238 :   for (j1 = j = 1; j1 < l1; j1++)
     248         189 :     if (gel(v1, j1)) gel(v,j++) = gel(v1,j1);
     249             : 
     250          49 :   LD = mkvecn(7, a1a2, b1b2, v, stoi(k), N, eno, r);
     251          49 :   if (!r) setlg(LD,7);
     252          49 :   return gerepilecopy(ltop, LD);
     253             : }
     254             : 
     255             : /*****************************************************************/
     256             : /*  L-series from closure                                        */
     257             : /*****************************************************************/
     258             : static GEN
     259       68474 : localfactor(void *E, GEN p, long n)
     260             : {
     261       68474 :   GEN s = closure_callgen2((GEN)E, p, utoi(n));
     262       68474 :   return direuler_factor(s, n);
     263             : }
     264             : static GEN
     265         665 : vecan_closure(GEN a, long L, long prec)
     266             : {
     267         665 :   long ta = typ(a);
     268         665 :   GEN gL, Sbad = NULL;
     269             : 
     270         665 :   if (ta == t_VEC)
     271             :   {
     272          70 :     long l = lg(a);
     273          70 :     if (l == 1) pari_err_TYPE("vecan_closure", a);
     274          70 :     ta = typ(gel(a,1));
     275             :     /* regular vector, return it */
     276          70 :     if (ta != t_CLOSURE) return vecslice(a, 1, minss(L,l-1));
     277          56 :     if (l != 3) pari_err_TYPE("vecan_closure", a);
     278          56 :     Sbad = gel(a,2);
     279          56 :     if (typ(Sbad) != t_VEC) pari_err_TYPE("vecan_closure", a);
     280          56 :     a = gel(a,1);
     281             :   }
     282         595 :   else if (ta != t_CLOSURE) pari_err_TYPE("vecan_closure", a);
     283         644 :   push_localprec(prec);
     284         644 :   gL = stoi(L);
     285         644 :   switch(closure_arity(a))
     286             :   {
     287             :     case 2:
     288         273 :       a = direuler_bad((void*)a, localfactor, gen_2, gL,gL, Sbad);
     289         273 :       break;
     290             :     case 1:
     291         371 :       a = closure_callgen1(a, gL);
     292         371 :       if (typ(a) != t_VEC) pari_err_TYPE("vecan_closure", a);
     293         371 :       break;
     294           0 :     default: pari_err_TYPE("vecan_closure [wrong arity]", a);
     295             :       a = NULL; /*LCOV_EXCL_LINE*/
     296             :   }
     297         644 :   pop_localprec(); return a;
     298             : }
     299             : 
     300             : /*****************************************************************/
     301             : /*  L-series of Dirichlet characters.                            */
     302             : /*****************************************************************/
     303             : 
     304             : static GEN
     305        1043 : lfunzeta(void)
     306             : {
     307        1043 :   GEN zet = mkvecn(7, NULL, gen_0, NULL, gen_1, gen_1, gen_1, gen_1);
     308        1043 :   gel(zet,1) = tag(gen_1, t_LFUN_ZETA);
     309        1043 :   gel(zet,3) = mkvec(gen_0);
     310        1043 :   return zet;
     311             : }
     312             : static GEN
     313         175 : lfunzetainit(GEN dom, long der, long bitprec)
     314         175 : { return lfuninit(lfunzeta(), dom, der, bitprec); }
     315             : 
     316             : static GEN
     317         574 : vecan_Kronecker(GEN D, long n)
     318             : {
     319         574 :   GEN v = cgetg(n+1, t_VECSMALL);
     320         574 :   ulong Du = itou_or_0(D);
     321         574 :   long i, id, d = Du ? minuu(Du, n): n;
     322         574 :   for (i = 1; i <= d; i++) v[i] = krois(D,i);
     323       19530 :   for (id = i; i <= n; i++,id++) /* periodic mod d */
     324             :   {
     325       18956 :     if (id > d) id = 1;
     326       18956 :     gel(v, i) = gel(v, id);
     327             :   }
     328         574 :   return v;
     329             : }
     330             : 
     331             : static GEN
     332         679 : lfunchiquad(GEN D)
     333             : {
     334             :   GEN r;
     335         679 :   if (equali1(D)) return lfunzeta();
     336         350 :   if (!isfundamental(D)) pari_err_TYPE("lfunchiquad [not primitive]", D);
     337         350 :   r = mkvecn(6, NULL, gen_0, NULL, gen_1, NULL, gen_1);
     338         350 :   gel(r,1) = tag(icopy(D), t_LFUN_KRONECKER);
     339         350 :   gel(r,3) = mkvec(signe(D) < 0? gen_1: gen_0);
     340         350 :   gel(r,5) = mpabs(D);
     341         350 :   return r;
     342             : }
     343             : 
     344             : /* Begin Hecke characters. Here a character is assumed to be given by a
     345             :    vector on the generators of the ray class group clgp of CL_m(K).
     346             :    If clgp = [h,[d1,...,dk],[g1,...,gk]] with dk|...|d2|d1, a character chi
     347             :    is given by [a1,a2,...,ak] such that chi(gi)=\zeta_di^ai. */
     348             : 
     349             : /* Value of CHI on x, coprime to bnr.mod */
     350             : static GEN
     351       28630 : chigeneval(GEN logx, GEN nchi, GEN z, long prec)
     352             : {
     353       28630 :   pari_sp av = avma;
     354       28630 :   GEN d = gel(nchi,1);
     355       28630 :   GEN e = FpV_dotproduct(gel(nchi,2), logx, d);
     356       28630 :   if (typ(z) != t_VEC)
     357           0 :     return gerepileupto(av, gpow(z, e, prec));
     358             :   else
     359             :   {
     360       28630 :     ulong i = itou(e);
     361       28630 :     avma = av; return gel(z, i+1);
     362             :   }
     363             : }
     364             : 
     365             : /* return x + yz; y != 0; z = 0,1 "often"; x = 0 "often" */
     366             : static GEN
     367      223538 : gaddmul(GEN x, GEN y, GEN z)
     368             : {
     369             :   pari_sp av;
     370      223538 :   if (typ(z) == t_INT)
     371             :   {
     372      211701 :     if (!signe(z)) return x;
     373       19278 :     if (equali1(z)) return gadd(x,y);
     374             :   }
     375       26068 :   if (isintzero(x)) return gmul(y,z);
     376        7812 :   av = avma;
     377        7812 :   return gerepileupto(av, gadd(x, gmul(y,z)));
     378             : }
     379             : 
     380             : static GEN
     381         294 : vecan_chiZ(GEN an, long n, long prec)
     382             : {
     383             :   forprime_t iter;
     384         294 :   GEN G = gel(an,1);
     385         294 :   GEN nchi = gel(an,2), gord = gel(nchi,1), z;
     386         294 :   GEN gp = cgetipos(3), v = vec_ei(n, 1);
     387         294 :   GEN N = znstar_get_N(G);
     388         294 :   long ord = itos_or_0(gord);
     389         294 :   ulong Nu = itou_or_0(N);
     390         294 :   long i, id, d = Nu ? minuu(Nu, n): n;
     391             :   ulong p;
     392             : 
     393         294 :   if (ord && n > (ord>>4))
     394         294 :   {
     395         294 :     GEN w = ncharvecexpo(G, nchi);
     396         294 :     z = grootsof1(ord, prec);
     397       10584 :     for (i = 1; i <= d; i++)
     398       10290 :       if (w[i] >= 0) gel(v, i) = gel(z, w[i]+1);
     399             :   }
     400             :   else
     401             :   {
     402           0 :     z = rootsof1_cx(gord, prec);
     403           0 :     u_forprime_init(&iter, 2, d);
     404           0 :     while ((p = u_forprime_next(&iter)))
     405             :     {
     406             :       GEN ch;
     407             :       ulong k;
     408           0 :       if (!umodiu(N,p)) continue;
     409           0 :       gp[2] = p;
     410           0 :       ch = chigeneval(znconreylog(G, gp), nchi, z, prec);
     411           0 :       gel(v, p)  = ch;
     412           0 :       for (k = 2*p; k <= (ulong)d; k += p)
     413           0 :         gel(v, k) = gaddmul(gel(v, k), ch, gel(v, k/p));
     414             :     }
     415             :   }
     416        2863 :   for (id = i = d+1; i <= n; i++,id++) /* periodic mod d */
     417             :   {
     418        2569 :     if (id > d) id = 1;
     419        2569 :     gel(v, i) = gel(v, id);
     420             :   }
     421         294 :   return v;
     422             : }
     423             : 
     424             : static GEN
     425         735 : vecan_chigen(GEN an, long n, long prec)
     426             : {
     427             :   forprime_t iter;
     428         735 :   GEN bnr = gel(an,1), nf = bnr_get_nf(bnr);
     429         735 :   GEN nchi = gel(an,2), gord = gel(nchi,1), z;
     430         735 :   GEN gp = cgetipos(3), v = vec_ei(n, 1);
     431         735 :   GEN N = gel(bnr_get_mod(bnr), 1), NZ = gcoeff(N,1,1);
     432         735 :   long ord = itos_or_0(gord);
     433             :   ulong p;
     434             : 
     435         735 :   if (ord && n > (ord>>4))
     436         735 :     z = grootsof1(ord, prec);
     437             :   else
     438           0 :     z = rootsof1_cx(gord, prec);
     439             : 
     440         735 :   if (nf_get_degree(nf) == 1)
     441             :   {
     442         553 :     ulong Nu = itou_or_0(NZ);
     443         553 :     long i, id, d = Nu ? minuu(Nu, n): n;
     444         553 :     u_forprime_init(&iter, 2, d);
     445        3276 :     while ((p = u_forprime_next(&iter)))
     446             :     {
     447             :       GEN ch;
     448             :       ulong k;
     449        2170 :       if (!umodiu(NZ,p)) continue;
     450        1610 :       gp[2] = p;
     451        1610 :       ch = chigeneval(isprincipalray(bnr,gp), nchi, z, prec);
     452        1610 :       gel(v, p)  = ch;
     453        3787 :       for (k = 2*p; k <= (ulong)d; k += p)
     454        2177 :         gel(v, k) = gaddmul(gel(v, k), ch, gel(v, k/p));
     455             :     }
     456        7490 :     for (id = i = d+1; i <= n; i++,id++) /* periodic mod d */
     457             :     {
     458        6937 :       if (id > d) id = 1;
     459        6937 :       gel(v, i) = gel(v, id);
     460             :     }
     461             :   }
     462             :   else
     463             :   {
     464         182 :     GEN BOUND = stoi(n);
     465         182 :     u_forprime_init(&iter, 2, n);
     466       27839 :     while ((p = u_forprime_next(&iter)))
     467             :     {
     468             :       GEN L;
     469             :       long j;
     470       27475 :       int check = !umodiu(NZ,p);
     471       27475 :       gp[2] = p;
     472       27475 :       L = idealprimedec_limit_norm(nf, gp, BOUND);
     473       54649 :       for (j = 1; j < lg(L); j++)
     474             :       {
     475       27174 :         GEN pr = gel(L, j), ch;
     476             :         ulong k, q;
     477       27174 :         if (check && idealval(nf, N, pr)) continue;
     478       27020 :         ch = chigeneval(isprincipalray(bnr,pr), nchi, z, prec);
     479       27020 :         q = upr_norm(pr);
     480       27020 :         gel(v, q) = gadd(gel(v, q), ch);
     481      248381 :         for (k = 2*q; k <= (ulong)n; k += q)
     482      221361 :           gel(v, k) = gaddmul(gel(v, k), ch, gel(v, k/q));
     483             :       }
     484             :     }
     485             :   }
     486         735 :   return v;
     487             : }
     488             : 
     489             : static GEN lfunzetak_i(GEN T);
     490             : static GEN
     491        3654 : vec01(long r1, long r2)
     492             : {
     493        3654 :   long d = r1+r2, i;
     494        3654 :   GEN v = cgetg(d+1,t_VEC);
     495        3654 :   for (i = 1; i <= r1; i++) gel(v,i) = gen_0;
     496        3654 :   for (     ; i <= d;  i++) gel(v,i) = gen_1;
     497        3654 :   return v;
     498             : }
     499             : 
     500             : /* G is a bid of nftyp typ_BIDZ */
     501             : static GEN
     502         210 : lfunchiZ(GEN G, GEN chi)
     503             : {
     504         210 :   pari_sp av = avma;
     505         210 :   GEN sig = NULL;
     506         210 :   GEN N = bid_get_ideal(G), nchi, r;
     507             :   int real;
     508             : 
     509         210 :   if (typ(N) != t_INT) pari_err_TYPE("lfunchiZ", G);
     510         210 :   if (equali1(N)) return lfunzeta();
     511         210 :   if (typ(chi) != t_COL) chi = znconreylog(G,chi);
     512         210 :   N = znconreyconductor(G, chi, &chi);
     513         210 :   if (typ(N) != t_INT)
     514             :   {
     515           0 :     if (equali1(gel(N,1))) { avma = av; return lfunzeta(); }
     516           0 :     G = znstar0(N, 1);
     517           0 :     N = gel(N,1);
     518             :   }
     519             :   /* chi now primitive on G */
     520         210 :   switch(itou_or_0(zncharorder(G, chi)))
     521             :   {
     522           0 :     case 1: avma = av; return lfunzeta();
     523          35 :     case 2: if (zncharisodd(G,chi)) N = negi(N);
     524          35 :             return gerepileupto(av, lfunchiquad(N));
     525             :   }
     526         175 :   sig = mkvec( zncharisodd(G, chi)? gen_1: gen_0 );
     527         175 :   nchi = znconreylog_normalize(G, chi);
     528         175 :   real = abscmpiu(gel(nchi,1), 2) <= 0;
     529         175 :   r = mkvecn(6, tag(mkvec2(G,nchi), t_LFUN_CHIZ),
     530             :                 real? gen_0: gen_1, sig, gen_1, N, gen_0);
     531         175 :   return gerepilecopy(av, r);
     532             : }
     533             : 
     534             : static GEN
     535        1148 : lfunchigen(GEN bnr, GEN CHI)
     536             : {
     537        1148 :   pari_sp av = avma;
     538             :   GEN v;
     539             :   GEN N, sig, Ldchi, nf, nchi, NN;
     540             :   long r1, r2, n1;
     541             :   int real;
     542             : 
     543        1148 :   if (nftyp(bnr) == typ_BIDZ) return lfunchiZ(bnr, CHI);
     544             : 
     545         945 :   v = bnrconductor_i(bnr, CHI, 2);
     546         945 :   bnr = gel(v,2);
     547         945 :   CHI = gel(v,3); /* now CHI is primitive wrt bnr */
     548             : 
     549         945 :   nf = bnr_get_nf(bnr);
     550         945 :   N = bnr_get_mod(bnr);
     551         945 :   n1 = lg(vec01_to_indices(gel(N,2))) - 1; /* vecsum(N[2]) */
     552         945 :   N = gel(N,1);
     553         945 :   NN = mulii(idealnorm(nf, N), absi(nf_get_disc(nf)));
     554         945 :   if (equali1(NN)) return gerepileupto(av, lfunzeta());
     555         602 :   if (ZV_equal0(CHI)) return gerepilecopy(av, lfunzetak_i(bnr));
     556         574 :   nf_get_sign(nf, &r1, &r2);
     557         574 :   sig = vec01(r1+r2-n1, r2+n1);
     558         574 :   nchi = char_normalize(CHI, cyc_normalize(bnr_get_cyc(bnr)));
     559         574 :   real = abscmpiu(gel(nchi,1), 2) <= 0;
     560         574 :   Ldchi = mkvecn(6, tag(mkvec2(bnr, nchi), t_LFUN_CHIGEN),
     561             :                     real? gen_0: gen_1, sig, gen_1, NN, gen_0);
     562         574 :   return gerepilecopy(av, Ldchi);
     563             : }
     564             : 
     565             : /* Find all characters of clgp whose kernel contain group given by HNF H.
     566             :  * Set *pcnj[i] if chi[i] is not real */
     567             : static GEN
     568         371 : chigenkerfind(GEN bnr, GEN H, GEN *pcnj)
     569             : {
     570         371 :   GEN res, cnj, L = bnrchar(bnr, H, NULL), cyc = bnr_get_cyc(bnr);
     571         371 :   long i, k, l = lg(L);
     572             : 
     573         371 :   res = cgetg(l, t_VEC);
     574         371 :   *pcnj = cnj = cgetg(l, t_VECSMALL);
     575        1414 :   for (i = k = 1; i < l; i++)
     576             :   {
     577        1043 :     GEN chi = gel(L,i), c = charconj(cyc, chi);
     578        1043 :     long fl = ZV_cmp(c, chi);
     579        1043 :     if (fl < 0) continue; /* keep one char in pair of conjugates */
     580         847 :     gel(res, k) = chi;
     581         847 :     cnj[k] = fl; k++;
     582             :   }
     583         371 :   setlg(cnj, k);
     584         371 :   setlg(res, k); return res;
     585             : }
     586             : 
     587             : static GEN
     588        1246 : lfunzetak_i(GEN T)
     589             : {
     590        1246 :   GEN Vga, N, nf, bnf = checkbnf_i(T), r = gen_0/*unknown*/;
     591             :   long r1, r2;
     592             : 
     593        1246 :   if (bnf)
     594         119 :     nf = bnf_get_nf(bnf);
     595             :   else
     596             :   {
     597        1127 :     nf = checknf_i(T);
     598        1127 :     if (!nf) nf = T = nfinit(T, DEFAULTPREC);
     599             :   }
     600        1246 :   nf_get_sign(nf,&r1,&r2);
     601        1246 :   N = absi(nf_get_disc(nf));
     602        1246 :   if (bnf)
     603             :   {
     604         119 :     GEN h = bnf_get_no(bnf);
     605         119 :     GEN R = bnf_get_reg(bnf);
     606         119 :     long prec = nf_get_prec(nf);
     607         238 :     r = gdiv(gmul(mulir(shifti(h, r1+r2), powru(mppi(prec), r2)), R),
     608         119 :              mulur(bnf_get_tuN(bnf), gsqrt(N, prec)));
     609             :   }
     610        1246 :   Vga = vec01(r1+r2,r2);
     611        1246 :   return mkvecn(7, tag(T,t_LFUN_NF), gen_0, Vga, gen_1, N, gen_1, r);
     612             : }
     613             : static GEN
     614         511 : lfunzetak(GEN T)
     615         511 : { pari_sp ltop = avma; return gerepilecopy(ltop, lfunzetak_i(T)); }
     616             : 
     617             : /* bnf = NULL: base field = Q */
     618             : GEN
     619         371 : lfunabelianrelinit(GEN nfabs, GEN bnf, GEN polrel, GEN dom, long der, long bitprec)
     620             : {
     621         371 :   pari_sp ltop = avma;
     622             :   GEN cond, chi, cnj, res, bnr, M, domain;
     623             :   long l, i;
     624         371 :   long v = -1;
     625             : 
     626         371 :   if (bnf) bnf = checkbnf(bnf);
     627             :   else
     628             :   {
     629         343 :     v = fetch_var();
     630         343 :     bnf = Buchall(pol_x(v), 0, nbits2prec(bitprec));
     631             :   }
     632         371 :   if (typ(polrel) != t_POL) pari_err_TYPE("lfunabelianrelinit", polrel);
     633         371 :   cond = rnfconductor(bnf, polrel);
     634         371 :   chi = chigenkerfind(gel(cond,2), gel(cond,3), &cnj);
     635         371 :   bnr = Buchray(bnf, gel(cond,1), nf_INIT);
     636         371 :   l = lg(chi); res = cgetg(l, t_VEC);
     637        1218 :   for (i = 1; i < l; ++i)
     638             :   {
     639         847 :     GEN L = lfunchigen(bnr, gel(chi,i));
     640         847 :     gel(res, i) = lfuninit(L, dom, der, bitprec);
     641             :   }
     642         371 :   if (v >= 0) delete_var();
     643         371 :   M = mkvec3(res, const_vecsmall(l-1, 1), cnj);
     644         371 :   domain = mkvec2(dom, mkvecsmall2(der, bitprec));
     645         371 :   return gerepilecopy(ltop, lfuninit_make(t_LDESC_PRODUCT, lfunzetak_i(nfabs), M, domain));
     646             : }
     647             : 
     648             : /*****************************************************************/
     649             : /*                 Dedekind zeta functions                       */
     650             : /*****************************************************************/
     651             : static GEN
     652        1239 : dirzetak0(GEN nf, ulong N)
     653             : {
     654        1239 :   GEN vect, c, c2, T = nf_get_pol(nf), index = nf_get_index(nf);
     655        1239 :   pari_sp av = avma, av2;
     656        1239 :   const ulong SQRTN = usqrt(N);
     657             :   ulong i, p, lx;
     658        1239 :   long court[] = {evaltyp(t_INT)|_evallg(3), evalsigne(1)|evallgefint(3),0};
     659             :   forprime_t S;
     660             : 
     661        1239 :   (void)evallg(N+1);
     662        1239 :   c  = cgetalloc(t_VECSMALL, N+1);
     663        1239 :   c2 = cgetalloc(t_VECSMALL, N+1);
     664        1239 :   c2[1] = c[1] = 1; for (i=2; i<=N; i++) c[i] = 0;
     665        1239 :   u_forprime_init(&S, 2, N);
     666        1239 :   av2 = avma;
     667      162239 :   while ( (p = u_forprime_next(&S)) )
     668             :   {
     669      159761 :     avma = av2;
     670      159761 :     if (umodiu(index, p)) /* p does not divide index */
     671      159558 :       vect = gel(Flx_degfact(ZX_to_Flx(T,p), p),1);
     672             :     else
     673             :     {
     674         203 :       court[2] = p;
     675         203 :       vect = idealprimedec_degrees(nf,court);
     676             :     }
     677      159761 :     lx = lg(vect);
     678      159761 :     if (p <= SQRTN)
     679       19341 :       for (i=1; i<lx; i++)
     680             :       {
     681       12866 :         ulong qn, q = upowuu(p, vect[i]); /* Norm P[i] */
     682       12866 :         if (!q || q > N) break;
     683       11011 :         memcpy(c2 + 2, c + 2, (N-1)*sizeof(long));
     684             :         /* c2[i] <- c[i] + sum_{k = 1}^{v_q(i)} c[i/q^k] for all i <= N */
     685       22750 :         for (qn = q; qn <= N; qn *= q)
     686             :         {
     687       22750 :           ulong k0 = N/qn, k, k2; /* k2 = k*qn */
     688       22750 :           for (k = k0, k2 = k*qn; k > 0; k--, k2 -=qn) c2[k2] += c[k];
     689       22750 :           if (q > k0) break; /* <=> q*qn > N */
     690             :         }
     691       11011 :         swap(c, c2);
     692             :       }
     693             :     else /* p > sqrt(N): simpler */
     694      298865 :       for (i=1; i<lx; i++)
     695             :       {
     696             :         ulong k, k2; /* k2 = k*p */
     697      261989 :         if (vect[i] > 1) break;
     698             :         /* c2[i] <- c[i] + sum_{k = 1}^{v_q(i)} c[i/q^k] for all i <= N */
     699      147434 :         for (k = N/p, k2 = k*p; k > 0; k--, k2 -= p) c[k2] += c[k];
     700             :       }
     701             :   }
     702        1239 :   avma = av;
     703        1239 :   pari_free(c2); return c;
     704             : }
     705             : 
     706             : GEN
     707        1239 : dirzetak(GEN nf, GEN b)
     708             : {
     709             :   GEN z, c;
     710             :   long n;
     711             : 
     712        1239 :   if (typ(b) != t_INT) pari_err_TYPE("dirzetak",b);
     713        1239 :   if (signe(b) <= 0) return cgetg(1,t_VEC);
     714        1239 :   nf = checknf(nf);
     715        1239 :   n = itou_or_0(b); if (!n) pari_err_OVERFLOW("dirzetak");
     716        1239 :   c = dirzetak0(nf, n);
     717        1239 :   z = vecsmall_to_vec(c); pari_free(c); return z;
     718             : }
     719             : 
     720             : static GEN
     721         672 : linit_get_mat(GEN linit)
     722             : {
     723         672 :   if (linit_get_type(linit)==t_LDESC_PRODUCT)
     724         161 :     return lfunprod_get_fact(linit_get_tech(linit));
     725             :   else
     726         511 :     return mkvec3(mkvec(linit), mkvecsmall(1), mkvecsmall(0));
     727             : }
     728             : 
     729             : static GEN
     730         336 : lfunproduct(GEN ldata, GEN linit1, GEN linit2, GEN domain)
     731             : {
     732         336 :   GEN M1 = linit_get_mat(linit1);
     733         336 :   GEN M2 = linit_get_mat(linit2);
     734        1344 :   GEN M3 = mkvec3(shallowconcat(gel(M1, 1), gel(M2, 1)),
     735         672 :                   vecsmall_concat(gel(M1, 2), gel(M2, 2)),
     736         672 :                   vecsmall_concat(gel(M1, 3), gel(M2, 3)));
     737         336 :   return lfuninit_make(t_LDESC_PRODUCT, ldata, M3, domain);
     738             : }
     739             : 
     740             : static GEN
     741           0 : lfunzetakinit_raw(GEN T, GEN dom, long der, long bitprec)
     742             : {
     743           0 :   pari_sp ltop = avma;
     744           0 :   GEN ldata = lfunzetak_i(T);
     745           0 :   return gerepileupto(ltop, lfuninit(ldata, dom, der, bitprec));
     746             : }
     747             : 
     748             : static GEN
     749         336 : lfunzetakinit_quotient(GEN nf, GEN polk, GEN dom, long der, long bitprec)
     750             : {
     751         336 :   pari_sp av = avma;
     752             :   GEN ak, an, nfk, Vga, ldata, N, Lk, LKk, domain;
     753             :   long r1k, r2k, r1, r2;
     754             : 
     755         336 :   nf_get_sign(nf,&r1,&r2);
     756         336 :   nfk = nfinit(polk, nbits2prec(bitprec));
     757         336 :   Lk = lfunzetakinit(nfk, dom, der, 0, bitprec); /* zeta_k */
     758         336 :   nf_get_sign(nfk,&r1k,&r2k);
     759         336 :   Vga = vec01((r1+r2) - (r1k+r2k), r2-r2k);
     760         336 :   N = absi(diviiexact(nf_get_disc(nf), nf_get_disc(nfk)));
     761         336 :   ak = nf_get_degree(nf)==1 ? tag(gen_1, t_LFUN_ZETA): tag(nfk, t_LFUN_NF);
     762         336 :   an = tag(mkvec2(tag(nf,t_LFUN_NF), ak), t_LFUN_DIV);
     763         336 :   ldata = mkvecn(6, an, gen_0, Vga, gen_1, N, gen_1);
     764         336 :   LKk = lfuninit(ldata, dom, der, bitprec); /* zeta_K/zeta_k */
     765         336 :   domain = mkvec2(dom, mkvecsmall2(der, bitprec));
     766         336 :   return gerepilecopy(av, lfunproduct(lfunzetak_i(nf), Lk, LKk, domain));
     767             : }
     768             : 
     769             : static GEN
     770          21 : subgroups_largestabelian(GEN S)
     771             : {
     772          21 :   long i, n = 0, l = lg (S);
     773          21 :   GEN M = NULL;
     774         147 :   for(i = 1; i < l; i++)
     775             :   {
     776         126 :     GEN Si = gel(S,i);
     777         126 :     long o = group_order(Si);
     778         126 :     if (o > n && group_isabelian(Si))
     779             :     {
     780          21 :       n = o;
     781          21 :       M = Si;
     782             :     }
     783             :   }
     784          21 :   return M;
     785             : }
     786             : 
     787             : 
     788             : /* If flag=0 (default), assume Artin conjecture */
     789             : 
     790             : static GEN
     791         364 : lfunzetakinit_Galois(GEN nf, GEN G, GEN dom, long der, long bitprec)
     792             : {
     793             :   GEN S, H, P, F, R, bnf;
     794         364 :   GEN T = nf_get_pol(nf);
     795         364 :   long v = varn(T);
     796         364 :   GEN grp = galois_group(G);
     797         364 :   if (group_isabelian(grp))
     798         343 :     return lfunabelianrelinit(nf, NULL, T, dom, der, bitprec);
     799          21 :   S = group_subgroups(grp);
     800          21 :   H = subgroups_largestabelian(S);
     801          21 :   if (v==0) { v=1; nf = gsubst(nf, 0, pol_x(v)); }
     802           7 :   else G = gsubst(G, v, pol_x(0));
     803          21 :   F = galoisfixedfield(G, H, 2, v);
     804          21 :   P = gel(F,1), R = gmael(F,3,1);
     805          21 :   setvarn(P, v);
     806          21 :   bnf = Buchall(P, 0, nbits2prec(bitprec));
     807          21 :   return lfunabelianrelinit(nf, bnf, R, dom, der, bitprec);
     808             : }
     809             : 
     810             : GEN
     811         875 : lfunzetakinit(GEN NF, GEN dom, long der, long flag, long bitprec)
     812             : {
     813         875 :   GEN nf = checknf(NF);
     814             :   GEN G, nfs, sbg;
     815         875 :   long lf, d = nf_get_degree(nf);
     816         875 :   if (d == 1) return lfunzetainit(dom, der, bitprec);
     817         700 :   G = galoisinit(nf, NULL);
     818         700 :   if (!isintzero(G))
     819         364 :     return lfunzetakinit_Galois(nf, G, dom, der, bitprec);
     820         336 :   nfs = nfsubfields(nf, 0); lf = lg(nfs)-1;
     821         336 :   sbg = gmael(nfs,lf-1,1);
     822         336 :   if (flag && d > 4*degpol(sbg))
     823           0 :     return lfunzetakinit_raw(nf, dom, der, bitprec);
     824         336 :   return lfunzetakinit_quotient(nf, sbg, dom, der, bitprec);
     825             : }
     826             : 
     827             : /***************************************************************/
     828             : /*             Elliptic Curves and Modular Forms               */
     829             : /***************************************************************/
     830             : 
     831             : static GEN
     832         168 : lfunellnf(GEN e)
     833             : {
     834         168 :   pari_sp av = avma;
     835             :   GEN ldata;
     836         168 :   GEN nf = ellnf_get_nf(e);
     837         168 :   GEN g = ellglobalred(e);
     838         168 :   GEN N = idealnorm(nf,gel(g,1)), d2 = sqri(nf_get_disc(nf));
     839         168 :   long n = nf_get_degree(nf);
     840         168 :   ldata = cgetg(7, t_VEC);
     841         168 :   gel(ldata, 1) = tag(e, t_LFUN_ELL);
     842         168 :   gel(ldata, 2) = gen_0;
     843         168 :   gel(ldata, 3) = vec01(n, n);
     844         168 :   gel(ldata, 4) = gen_2;
     845         168 :   gel(ldata, 5) = mulii(d2,N);
     846         168 :   gel(ldata, 6) = stoi(ellrootno_global(e));
     847         168 :   return gerepileupto(av, ldata);
     848             : }
     849             : 
     850             : static GEN
     851         609 : lfunellQ(GEN e)
     852             : {
     853         609 :   pari_sp av = avma;
     854         609 :   GEN ldata = cgetg(7, t_VEC);
     855         609 :   gel(ldata, 1) = tag(ellanal_globalred(e, NULL), t_LFUN_ELL);
     856         609 :   gel(ldata, 2) = gen_0;
     857         609 :   gel(ldata, 3) = mkvec2(gen_0, gen_1);
     858         609 :   gel(ldata, 4) = gen_2;
     859         609 :   gel(ldata, 5) = icopy(ellQ_get_N(e));
     860         609 :   gel(ldata, 6) = stoi(ellrootno_global(e));
     861         609 :   return gerepilecopy(av, ldata); /* ellanal_globalred not gerepile-safe */
     862             : }
     863             : 
     864             : static GEN
     865         777 : lfunell(GEN e)
     866             : {
     867         777 :   long t = ell_get_type(e);
     868         777 :   switch(t)
     869             :   {
     870             :     case t_ELL_Q:
     871         609 :       return lfunellQ(e);
     872             :     case t_ELL_NF:
     873         168 :       return lfunellnf(e);
     874             :   }
     875           0 :   pari_err_TYPE("lfun",e);
     876             :   return NULL; /*LCOV_EXCL_LINE*/
     877             : }
     878             : 
     879             : GEN
     880           7 : lfunmfspec(GEN lmisc, long bitprec)
     881             : {
     882           7 :   pari_sp ltop = avma;
     883             :   GEN Vga, linit, ldataf, veven, vodd, om, op, eps, dom;
     884             :   long k, k2, j;
     885             : 
     886           7 :   ldataf = lfunmisc_to_ldata_shallow(lmisc);
     887           7 :   k = ldata_get_k(ldataf);
     888           7 :   dom = mkvec3(dbltor(k/2.), dbltor((k-2)/2.), gen_0);
     889           7 :   if (is_linit(lmisc) && linit_get_type(lmisc) == t_LDESC_INIT
     890           0 :       && sdomain_isincl(k, dom, lfun_get_dom(linit_get_tech(lmisc))))
     891           0 :     linit = lmisc;
     892             :   else
     893           7 :     linit = lfuninit(ldataf, dom, 0, bitprec);
     894           7 :   Vga = ldata_get_gammavec(ldataf);
     895           7 :   if (!ldata_isreal(ldataf) || !gequal(Vga, mkvec2(gen_0,gen_1)))
     896           0 :     pari_err_TYPE("lfunmfspec", lmisc);
     897           7 :   if (odd(k)) pari_err_IMPL("odd weight in lfunmfspec");
     898           7 :   k2 = k/2;
     899           7 :   vodd = cgetg(k2+1, t_VEC);
     900           7 :   veven = cgetg(k2, t_VEC);
     901           7 :   for (j=1; j <= k2; ++j) gel(vodd,j) = lfunlambda(linit, stoi(2*j-1), bitprec);
     902           7 :   for (j=1; j < k2; ++j) gel(veven,j) = lfunlambda(linit, stoi(2*j), bitprec);
     903           7 :   if (k > 2)
     904             :   {
     905           7 :     om = gel(veven,1);
     906           7 :     veven = gdiv(veven, om);
     907           7 :     op = gel(vodd,2);
     908             :   }
     909             :   else
     910             :   { /* veven empty */
     911           0 :     om = gen_1;
     912           0 :     op = gel(vodd,1);
     913             :   }
     914           7 :   vodd = gdiv(vodd, op);
     915           7 :   eps = int2n(bitprec/4);
     916           7 :   veven= bestappr(veven, eps);
     917           7 :   vodd = bestappr(vodd, eps);
     918           7 :   return gerepilecopy(ltop, mkvec4(veven, vodd, om, op));
     919             : }
     920             : 
     921             : static long
     922          42 : ellsymsq_bad2(GEN c4, GEN c6, long e, long *pb2)
     923             : {
     924          42 :   switch (e)
     925             :   {
     926          14 :     case 2: *pb2 = 1; return 1;
     927          14 :     case 3: *pb2 = 2; return 0;
     928          14 :     case 5: *pb2 = 3; return 0;
     929           0 :     case 7: *pb2 = 4; return 0;
     930             :     case 8:
     931           0 :       if (dvdiu(c6,512)) { *pb2 = 4; return 0; }
     932           0 :       *pb2 = 3; return umodiu(c4,128)==32 ? 1 : -1;
     933           0 :     default: *pb2 = 0; return 0;
     934             :   }
     935             : }
     936             : static long
     937          28 : ellsymsq_bad3(GEN c4, GEN c6, long e, long *pb3)
     938             : {
     939             :   long c6_243, c4_81;
     940          28 :   switch (e)
     941             :   {
     942           0 :     case 2: *pb3 = 1; return 1;
     943          28 :     case 3: *pb3 = 2; return 0;
     944           0 :     case 5: *pb3 = 3; return 0;
     945           0 :     case 4: *pb3 = 2;
     946           0 :       c4_81 = umodiu(c4,81);
     947           0 :       if (c4_81 == 27) return -1;
     948           0 :       if (c4_81%27 != 9) return 1;
     949           0 :       c6_243 = umodiu(c6,243);
     950           0 :       return (c6_243==108 || c6_243==135)? -1: 1;
     951           0 :     default: *pb3 = 0; return 0;
     952             :   }
     953             : }
     954             : static int
     955           0 : c4c6_testp(GEN c4, GEN c6, GEN p)
     956           0 : { GEN p2 = sqri(p); return (dvdii(c6,p2) && !dvdii(c4,p2)); }
     957             : /* assume e = v_p(N) >= 2 */
     958             : static long
     959          70 : ellsymsq_badp(GEN c4, GEN c6, GEN p, long e, long *pb)
     960             : {
     961          70 :   if (absequaliu(p, 2)) return ellsymsq_bad2(c4, c6, e, pb);
     962          28 :   if (absequaliu(p, 3)) return ellsymsq_bad3(c4, c6, e, pb);
     963           0 :   *pb = 1;
     964           0 :   switch(umodiu(p, 12UL))
     965             :   {
     966           0 :     case 1: return -1;
     967           0 :     case 5: return c4c6_testp(c4,c6,p)? -1: 1;
     968           0 :     case 7: return c4c6_testp(c4,c6,p)?  1:-1;
     969           0 :     default:return 1; /* p%12 = 11 */
     970             :   }
     971             : }
     972             : static GEN
     973        5012 : ellsymsq(void *D, GEN p, long n)
     974             : {
     975        5012 :   GEN E = (GEN)D;
     976        5012 :   GEN T, ap = sqri(ellap(E, p));
     977        5012 :   long e = Z_pval(ellQ_get_N(E), p);
     978        5012 :   if (e)
     979             :   {
     980          91 :     if (e == 1)
     981          56 :       T = deg1pol_shallow(negi(ap),gen_1,0);
     982             :     else
     983             :     {
     984          35 :       GEN c4 = ell_get_c4(E);
     985          35 :       GEN c6 = ell_get_c6(E);
     986          35 :       long junk, a = ellsymsq_badp(c4, c6, p, e, &junk);
     987          35 :       GEN pb = negi(mulis(p,a));
     988          35 :       GEN u1 = negi(addii(ap,pb));
     989          35 :       GEN u2 = mulii(ap,pb);
     990          35 :       T = mkpoln(3,u2,u1,gen_1);
     991             :     }
     992             :   }
     993             :   else
     994             :   {
     995        4921 :     GEN u1 = subii(ap,p);
     996        4921 :     GEN u2 = mulii(p,u1);
     997        4921 :     GEN u3 = powiu(p,3);
     998        4921 :     T = mkpoln(4,negi(u3),u2,negi(u1),gen_1);
     999             :   }
    1000        5012 :   return RgXn_inv(T,n);
    1001             : }
    1002             : static GEN
    1003          56 : vecan_ellsymsq(GEN an, long n)
    1004          56 : { GEN nn = stoi(n); return direuler_bad((void*)an, &ellsymsq, gen_2, nn, nn, NULL); }
    1005             : 
    1006             : static GEN
    1007          56 : lfunellsymsqmintwist(GEN e)
    1008             : {
    1009          56 :   pari_sp av = avma;
    1010             :   GEN B, N, Nfa, P, E, V, c4, c6, ld;
    1011             :   long i, l, k;
    1012          56 :   checkell_Q(e);
    1013          56 :   e = ellminimalmodel(e, NULL);
    1014          56 :   ellQ_get_Nfa(e, &N, &Nfa);
    1015          56 :   c4 = ell_get_c4(e);
    1016          56 :   c6 = ell_get_c6(e);
    1017          56 :   P = gel(Nfa,1); l = lg(P);
    1018          56 :   E = gel(Nfa,2);
    1019          56 :   V = cgetg(l, t_VEC);
    1020          56 :   B = gen_1;
    1021         147 :   for (i=k=1; i<l; i++)
    1022             :   {
    1023          91 :     GEN p = gel(P,i);
    1024          91 :     long a, b, e = itos(gel(E,i));
    1025          91 :     if (e == 1) { B = mulii(B, p); continue; }
    1026          35 :     a = ellsymsq_badp(c4, c6, p, e, &b);
    1027          35 :     B = mulii(B, powiu(p, b));
    1028          35 :     gel(V,k++) = mkvec2(p, stoi(a));
    1029             :   }
    1030          56 :   setlg(V, k);
    1031          56 :   ld = mkvecn(6, tag(e, t_LFUN_SYMSQ_ELL), gen_0,
    1032             :                  mkvec3(gen_0, gen_0, gen_1), stoi(3), sqri(B), gen_1);
    1033          56 :   return gerepilecopy(av, mkvec2(ld, V));
    1034             : }
    1035             : 
    1036             : static GEN
    1037          56 : mfpeters(GEN ldata2, GEN fudge, GEN N, long k, long bitprec)
    1038             : {
    1039          56 :   GEN t, L = greal(lfun(ldata2, stoi(k), bitprec));
    1040          56 :   long prec = nbits2prec(bitprec);
    1041          56 :   t = powrs(mppi(prec), k+1); shiftr_inplace(t, 2*k-1); /* Pi/2 * (4Pi)^k */
    1042          56 :   return gmul(gdiv(gmul(mulii(N,mpfact(k-1)), fudge), t), L);
    1043             : }
    1044             : 
    1045             : /* Assume E to be twist-minimal */
    1046             : static GEN
    1047          56 : lfunellmfpetersmintwist(GEN E, long bitprec)
    1048             : {
    1049          56 :   pari_sp av = avma;
    1050          56 :   GEN symsq, veceuler, N = ellQ_get_N(E), fudge = gen_1;
    1051          56 :   long j, k = 2;
    1052          56 :   symsq = lfunellsymsqmintwist(E);
    1053          56 :   veceuler = gel(symsq,2);
    1054          91 :   for (j = 1; j < lg(veceuler); j++)
    1055             :   {
    1056          35 :     GEN v = gel(veceuler,j), p = gel(v,1), q = powis(p,1-k);
    1057          35 :     long s = signe(gel(v,2));
    1058          35 :     if (s) fudge = gmul(fudge, s==1 ? gaddsg(1, q): gsubsg(1, q));
    1059             :   }
    1060          56 :   return gerepileupto(av, mfpeters(gel(symsq,1),fudge,N,k,bitprec));
    1061             : }
    1062             : 
    1063             : /* From Christophe Delaunay, http://delaunay.perso.math.cnrs.fr/these.pdf */
    1064             : static GEN
    1065          56 : elldiscfix(GEN E, GEN Et, GEN D)
    1066             : {
    1067          56 :   GEN N = ellQ_get_N(E), Nt = ellQ_get_N(Et);
    1068          56 :   GEN P = gel(Z_factor(absi(D)), 1);
    1069          56 :   GEN f = gen_1;
    1070          56 :   long i, l = lg(P);
    1071         105 :   for (i=1; i < l; i++)
    1072             :   {
    1073          49 :     GEN r, p = gel(P,i);
    1074          49 :     long v = Z_pval(N, p), vt = Z_pval(Nt, p);
    1075          49 :     if (v <= vt) continue;
    1076             :     /* v > vt */
    1077          35 :     if (absequaliu(p, 2))
    1078             :     {
    1079          28 :       if (vt == 0 && v >= 4)
    1080           0 :         r = shifti(subsi(9, sqri(ellap(Et, p))), v-3);  /* 9=(2+1)^2 */
    1081          28 :       else if (vt == 1)
    1082           7 :         r = gmul2n(utoipos(3), v-3);  /* not in Z if v=2 */
    1083          21 :       else if (vt >= 2)
    1084          21 :         r = int2n(v-vt);
    1085             :       else
    1086           0 :         r = gen_1; /* vt = 0, 1 <= v <= 3 */
    1087             :     }
    1088           7 :     else if (vt >= 1)
    1089           7 :       r = gdiv(subiu(sqri(p), 1), p);
    1090             :     else
    1091           0 :       r = gdiv(mulii(subiu(p, 1), subii(sqri(addiu(p, 1)), sqri(ellap(Et, p)))), p);
    1092          35 :     f = gmul(f, r);
    1093             :   }
    1094          56 :   return f;
    1095             : }
    1096             : 
    1097             : GEN
    1098          56 : lfunellmfpeters(GEN E, long bitprec)
    1099             : {
    1100          56 :   pari_sp ltop = avma;
    1101          56 :   long prec = nbits2prec(bitprec);
    1102          56 :   GEN D = ellminimaltwistcond(E);
    1103          56 :   GEN Etr = ellinit(elltwist(E, D), NULL, prec);
    1104          56 :   GEN Et = ellminimalmodel(Etr, NULL);
    1105          56 :   GEN nor = lfunellmfpetersmintwist(Et, bitprec);
    1106          56 :   GEN nor2 = gmul(nor, elldiscfix(E, Et, D));
    1107          56 :   obj_free(Etr); obj_free(Et);
    1108          56 :   return gerepilecopy(ltop, nor2);
    1109             : }
    1110             : 
    1111             : /*************************************************************/
    1112             : /*               Genus 2 curves                              */
    1113             : /*************************************************************/
    1114             : 
    1115             : static void
    1116      187446 : Flv_diffnext(GEN d, ulong p)
    1117             : {
    1118      187446 :   long j, n = lg(d)-1;
    1119     1312122 :   for(j = n; j>=2; j--)
    1120     1124676 :     uel(d,j) = Fl_add(uel(d,j), uel(d,j-1), p);
    1121      187446 : }
    1122             : 
    1123             : static GEN
    1124        1386 : Flx_difftable(GEN P, ulong p)
    1125             : {
    1126        1386 :   long i, n = degpol(P);
    1127        1386 :   GEN V = cgetg(n+2, t_VECSMALL);
    1128        1386 :   uel(V, n+1) = uel(P, 2);
    1129        9702 :   for(i = n; i >= 1; i--)
    1130             :   {
    1131        8316 :     P = Flx_diff1(P, p);
    1132        8316 :     uel(V, i) = uel(P, 2);
    1133             :   }
    1134        1386 :   return V;
    1135             : }
    1136             : 
    1137             : static long
    1138        1386 : Flx_genus2trace_naive(GEN H, ulong p)
    1139             : {
    1140        1386 :   pari_sp av = avma;
    1141             :   ulong i, j;
    1142        1386 :   long a, n = degpol(H);
    1143        1386 :   GEN k = const_vecsmall(p, -1), d;
    1144        1386 :   k[1] = 0;
    1145       94416 :   for (i=1, j=1; i < p; i += 2, j = Fl_add(j, i, p))
    1146       93030 :     k[j+1] = 1;
    1147        1386 :   a = n == 5 ? 0: k[1+Flx_lead(H)];
    1148        1386 :   d = Flx_difftable(H, p);
    1149      188832 :   for (i=0; i < p; i++)
    1150             :   {
    1151      187446 :     a += k[1+uel(d,n+1)];
    1152      187446 :     Flv_diffnext(d, p);
    1153             :   }
    1154        1386 :   avma = av;
    1155        1386 :   return a;
    1156             : }
    1157             : 
    1158             : static GEN
    1159        1512 : dirgenus2(void *E, GEN p, long n)
    1160             : {
    1161        1512 :   pari_sp av = avma;
    1162        1512 :   GEN Q = (GEN) E;
    1163             :   GEN f;
    1164        1512 :   if (n > 2)
    1165         126 :     f = RgX_recip(hyperellcharpoly(gmul(Q,gmodulo(gen_1, p))));
    1166             :   else
    1167             :   {
    1168        1386 :     ulong pp = itou(p);
    1169        1386 :     GEN Qp = ZX_to_Flx(Q, pp);
    1170        1386 :     long t = Flx_genus2trace_naive(Qp, pp);
    1171        1386 :     f = deg1pol_shallow(stoi(t), gen_1, 0);
    1172             :   }
    1173        1512 :   return gerepileupto(av, RgXn_inv(f, n));
    1174             : }
    1175             : 
    1176             : static GEN
    1177          28 : vecan_genus2(GEN an, long L)
    1178             : {
    1179          28 :   GEN Q = gel(an,1), bad = gel(an, 2);
    1180          28 :   return direuler_bad((void*)Q, dirgenus2, gen_2, stoi(L), NULL, bad);
    1181             : }
    1182             : 
    1183             : static GEN
    1184          14 : genus2_redmodel(GEN P, GEN p)
    1185             : {
    1186          14 :   GEN M = FpX_factor(P, p);
    1187          14 :   GEN F = gel(M,1), E = gel(M,2);
    1188          14 :   long i, k, r = lg(F);
    1189          14 :   GEN U = scalarpol(leading_coeff(P), varn(P));
    1190          14 :   GEN G = cgetg(r, t_COL);
    1191          49 :   for (i=1, k=0; i<r; i++)
    1192             :   {
    1193          35 :     if (E[i]>1)
    1194          28 :       gel(G,++k) = gel(F,i);
    1195          35 :     if (odd(E[i]))
    1196          14 :       U = FpX_mul(U, gel(F,i), p);
    1197             :   }
    1198          14 :   setlg(G,++k);
    1199          14 :   return mkvec2(G,U);
    1200             : }
    1201             : 
    1202             : static GEN
    1203         140 : oneminusxd(long d)
    1204             : {
    1205         140 :   return gsub(gen_1, pol_xn(d, 0));
    1206             : }
    1207             : 
    1208             : static GEN
    1209           7 : ellfromeqncharpoly(GEN P, GEN Q, GEN p)
    1210             : {
    1211             :   long v;
    1212             :   GEN E, F, t, y;
    1213           7 :   v = fetch_var();
    1214           7 :   y = pol_x(v);
    1215           7 :   F = gsub(gadd(ZX_sqr(y), gmul(y, Q)), P);
    1216           7 :   E = ellinit(ellfromeqn(F), p, DEFAULTPREC);
    1217           7 :   delete_var();
    1218           7 :   t = ellap(E, p);
    1219           7 :   obj_free(E);
    1220           7 :   return mkpoln(3, gen_1, negi(t), p);
    1221             : }
    1222             : 
    1223             : static GEN
    1224          14 : genus2_eulerfact(GEN P, GEN p)
    1225             : {
    1226          14 :   GEN Pp = FpX_red(P, p);
    1227          14 :   GEN GU = genus2_redmodel(Pp, p);
    1228          14 :   long d = 6-degpol(Pp), v = d/2, w = odd(d);
    1229             :   GEN abe, tor;
    1230          14 :   GEN ki, kp = pol_1(0), kq = pol_1(0);
    1231          14 :   GEN F = gel(GU,1), Q = gel(GU,2);
    1232          14 :   long dQ = degpol(Q), lF = lg(F)-1;
    1233             : 
    1234          14 :   abe = dQ >= 5 ? RgX_recip(hyperellcharpoly(gmul(Q,gmodulo(gen_1,p))))
    1235          35 :       : dQ >= 3 ? RgX_recip(ellfromeqncharpoly(Q,gen_0,p))
    1236          21 :                 : pol_1(0);
    1237          14 :   ki = dQ > 0 ? oneminusxd(1)
    1238          28 :               : Fp_issquare(gel(Q,2),p) ? ZX_sqr(oneminusxd(1))
    1239          14 :                                         : oneminusxd(2);
    1240          14 :   if (lF)
    1241             :   {
    1242             :     long i;
    1243          42 :     for(i=1; i <= lF; i++)
    1244             :     {
    1245          28 :       GEN Fi = gel(F, i);
    1246          28 :       long d = degpol(Fi);
    1247          28 :       GEN e = FpX_rem(Q, Fi, p);
    1248          49 :       GEN kqf = lgpol(e)==0 ? oneminusxd(d):
    1249          42 :                 FpXQ_issquare(e, Fi, p) ? ZX_sqr(oneminusxd(d))
    1250          42 :                                         : oneminusxd(2*d);
    1251          28 :       kp = gmul(kp, oneminusxd(d));
    1252          28 :       kq = gmul(kq, kqf);
    1253             :     }
    1254             :   }
    1255          14 :   if (v)
    1256             :   {
    1257           0 :     GEN kqoo = w==1 ? oneminusxd(1):
    1258           0 :                Fp_issquare(leading_coeff(Q), p)? ZX_sqr(oneminusxd(1))
    1259           0 :                                               : oneminusxd(2);
    1260           0 :     kp = gmul(kp, oneminusxd(1));
    1261           0 :     kq = gmul(kq, kqoo);
    1262             :   }
    1263          14 :   tor = RgX_div(ZX_mul(oneminusxd(1), kq), ZX_mul(ki, kp));
    1264          14 :   return ginv( ZX_mul(abe, tor) );
    1265             : }
    1266             : 
    1267             : static GEN
    1268          14 : F2x_genus2_find_trans(GEN P, GEN Q, GEN F)
    1269             : {
    1270          14 :   pari_sp av = avma;
    1271          14 :   long i, d = F2x_degree(F), v = P[1];
    1272             :   GEN M, C, V;
    1273          14 :   M = cgetg(d+1, t_MAT);
    1274          42 :   for (i=1; i<=d; i++)
    1275             :   {
    1276          28 :     GEN Mi = F2x_rem(F2x_add(F2x_shift(Q,i-1), monomial_F2x(2*i-2,v)), F);
    1277          28 :     gel(M,i) = F2x_to_F2v(Mi, d);
    1278             :   }
    1279          14 :   C = F2x_to_F2v(F2x_rem(P, F), d);
    1280          14 :   V = F2m_F2c_invimage(M, C);
    1281          14 :   return gerepileuptoleaf(av, F2v_to_F2x(V, v));
    1282             : }
    1283             : 
    1284             : static GEN
    1285          21 : F2x_genus2_trans(GEN P, GEN Q, GEN H)
    1286             : {
    1287          21 :   return F2x_add(P,F2x_add(F2x_mul(H,Q), F2x_sqr(H)));
    1288             : }
    1289             : 
    1290             : static GEN
    1291          21 : F2x_genus_redoo(GEN P, GEN Q, long k)
    1292             : {
    1293          21 :   if (F2x_degree(P)==2*k)
    1294             :   {
    1295          14 :     long c = F2x_coeff(P,2*k-1), dQ = F2x_degree(Q);
    1296          14 :     if ((dQ==k-1 && c==1) || (dQ<k-1 && c==0))
    1297           7 :      return F2x_genus2_trans(P, Q, monomial_F2x(k, P[1]));
    1298             :   }
    1299          14 :   return P;
    1300             : }
    1301             : 
    1302             : static GEN
    1303          14 : F2x_pseudodisc(GEN P, GEN Q)
    1304             : {
    1305          14 :   GEN dP = F2x_deriv(P), dQ = F2x_deriv(Q);
    1306          14 :   return F2x_gcd(Q, F2x_add(F2x_mul(P, F2x_sqr(dQ)), F2x_sqr(dP)));
    1307             : }
    1308             : 
    1309             : static GEN
    1310           7 : F2x_genus_red(GEN P, GEN Q)
    1311             : {
    1312             :   long dP, dQ;
    1313             :   GEN F, FF;
    1314           7 :   P = F2x_genus_redoo(P, Q, 3);
    1315           7 :   P = F2x_genus_redoo(P, Q, 2);
    1316           7 :   P = F2x_genus_redoo(P, Q, 1);
    1317           7 :   dP = F2x_degree(P);
    1318           7 :   dQ = F2x_degree(Q);
    1319           7 :   FF = F = F2x_pseudodisc(P,Q);
    1320          21 :   while(F2x_degree(F)>0)
    1321             :   {
    1322           7 :     GEN M = gel(F2x_factor(F),1);
    1323           7 :     long i, l = lg(M);
    1324          21 :     for(i=1; i<l; i++)
    1325             :     {
    1326          14 :       GEN R = F2x_sqr(gel(M,i));
    1327          14 :       GEN H = F2x_genus2_find_trans(P, Q, R);
    1328          14 :       P = F2x_div(F2x_genus2_trans(P, Q, H), R);
    1329          14 :       Q = F2x_div(Q, gel(M,i));
    1330             :     }
    1331           7 :     F = F2x_pseudodisc(P, Q);
    1332             :   }
    1333           7 :   return mkvec4(P,Q,FF,mkvecsmall2(dP,dQ));
    1334             : }
    1335             : 
    1336             : /* Number of solutions of x^2+b*x+c */
    1337             : static long
    1338          14 : F2xqX_quad_nbroots(GEN b, GEN c, GEN T)
    1339             : {
    1340          14 :   if (lgpol(b) > 0)
    1341             :   {
    1342          14 :     GEN d = F2xq_div(c, F2xq_sqr(b, T), T);
    1343          14 :     return F2xq_trace(d, T)? 0: 2;
    1344             :   }
    1345             :   else
    1346           0 :     return 1;
    1347             : }
    1348             : 
    1349             : static GEN
    1350          14 : genus2_redmodel2(GEN P)
    1351             : {
    1352          14 :   GEN Q = pol_0(varn(P));
    1353          14 :   GEN P2 = ZX_to_F2x(P);
    1354          42 :   while (F2x_issquare(P2))
    1355             :   {
    1356          14 :     GEN H = F2x_to_ZX(F2x_sqrt(P2));
    1357          14 :     GEN P1 = ZX_sub(P, ZX_sqr(H));
    1358          14 :     GEN Q1 = ZX_add(Q, ZX_mulu(H, 2));
    1359          14 :     if ((signe(P1)==0 ||  ZX_lval(P1,2)>=2)
    1360          14 :      && (signe(Q1)==0 ||  ZX_lval(Q1,2)>=1))
    1361             :     {
    1362          14 :       P = ZX_shifti(P1, -2);
    1363          14 :       Q = ZX_shifti(Q1, -1);
    1364          14 :       P2= ZX_to_F2x(P);
    1365             :     } else break;
    1366             :   }
    1367          14 :   return mkvec2(P,Q);
    1368             : }
    1369             : 
    1370             : static GEN
    1371           7 : genus2_eulerfact2(GEN PQ)
    1372             : {
    1373           7 :   GEN V = F2x_genus_red(ZX_to_F2x(gel(PQ, 1)), ZX_to_F2x(gel(PQ, 2)));
    1374           7 :   GEN P = gel(V, 1), Q = gel(V, 2);
    1375           7 :   GEN F = gel(V, 3), v = gel(V, 4);
    1376             :   GEN abe, tor;
    1377           7 :   GEN ki, kp = pol_1(0), kq = pol_1(0);
    1378           7 :   long dP = F2x_degree(P), dQ = F2x_degree(Q), d = maxss(dP, 2*dQ);
    1379           7 :   if (!lgpol(F)) return pol_1(0);
    1380          14 :   ki = dQ!=0 || dP>0 ? oneminusxd(1):
    1381           7 :       dP==-1 ? ZX_sqr(oneminusxd(1)): oneminusxd(2);
    1382          14 :   abe = d>=5? RgX_recip(hyperellcharpoly(gmul(PQ,gmodulss(1,2)))):
    1383           7 :         d>=3? RgX_recip(ellfromeqncharpoly(F2x_to_ZX(P), F2x_to_ZX(Q), gen_2)):
    1384             :         pol_1(0);
    1385           7 :   if (lgpol(F))
    1386             :   {
    1387           7 :     GEN M = gel(F2x_factor(F), 1);
    1388           7 :     long i, lF = lg(M)-1;
    1389          21 :     for(i=1; i <= lF; i++)
    1390             :     {
    1391          14 :       GEN Fi = gel(M, i);
    1392          14 :       long d = F2x_degree(Fi);
    1393          14 :       long nb  = F2xqX_quad_nbroots(F2x_rem(Q, Fi), F2x_rem(P, Fi), Fi);
    1394          28 :       GEN kqf = nb==1 ? oneminusxd(d):
    1395           0 :                 nb==2 ? ZX_sqr(oneminusxd(d))
    1396          14 :                       : oneminusxd(2*d);
    1397          14 :       kp = gmul(kp, oneminusxd(d));
    1398          14 :       kq = gmul(kq, kqf);
    1399             :     }
    1400             :   }
    1401           7 :   if (maxss(v[1],2*v[2])<5)
    1402             :   {
    1403          14 :     GEN kqoo = v[1]>2*v[2] ? oneminusxd(1):
    1404           7 :                v[1]<2*v[2] ? ZX_sqr(oneminusxd(1))
    1405           7 :                            : oneminusxd(2);
    1406           7 :     kp = gmul(kp, oneminusxd(1));
    1407           7 :     kq = gmul(kq, kqoo);
    1408             :   }
    1409           7 :   tor = RgX_div(ZX_mul(oneminusxd(1),kq), ZX_mul(ki, kp));
    1410           7 :   return ginv( ZX_mul(abe, tor) );
    1411             : }
    1412             : 
    1413             : GEN
    1414          14 : lfungenus2(GEN G)
    1415             : {
    1416          14 :   pari_sp ltop = avma;
    1417             :   GEN Ldata;
    1418          14 :   GEN gr = genus2red(G, NULL);
    1419          14 :   GEN N  = gel(gr, 1), M = gel(gr, 2), Q = gel(gr, 3), L = gel(gr, 4);
    1420          14 :   GEN PQ = genus2_redmodel2(Q);
    1421             :   GEN e;
    1422          14 :   long i, lL = lg(L), ram2;
    1423          14 :   ram2 = absequaliu(gmael(M,1,1),2);
    1424          14 :   if (ram2 && equalis(gmael(M,2,1),-1))
    1425           7 :     pari_warn(warner,"unknown valuation of conductor at 2");
    1426          14 :   e = cgetg(lL+(ram2?0:1), t_VEC);
    1427          21 :   gel(e,1) = mkvec2(gen_2, ram2 ? genus2_eulerfact2(PQ)
    1428           7 :            : ginv( RgX_recip(hyperellcharpoly(gmul(PQ,gmodulss(1,2))))) );
    1429          28 :   for(i = ram2? 2: 1; i < lL; i++)
    1430             :   {
    1431          14 :     GEN Li = gel(L, i);
    1432          14 :     GEN p = gel(Li, 1);
    1433          14 :     gel(e, ram2 ? i: i+1) = mkvec2(p, genus2_eulerfact(Q,p));
    1434             :   }
    1435          14 :   Ldata = mkvecn(6, tag(mkvec2(Q,e), t_LFUN_GENUS2),
    1436             :       gen_0, mkvec4(gen_0, gen_0, gen_1, gen_1), gen_2, N, gen_0);
    1437          14 :   return gerepilecopy(ltop, Ldata);
    1438             : }
    1439             : 
    1440             : /*************************************************************/
    1441             : /*                        ETA QUOTIENTS                      */
    1442             : /* An eta quotient is a matrix with 2 columns [m, r_m] with  */
    1443             : /* m >= 1 representing f(\tau)=\prod_m\eta(m\tau)^{r_m}.     */
    1444             : /*************************************************************/
    1445             : 
    1446             : GEN
    1447         476 : eta_inflate_ZXn(long m, long v)
    1448             : {
    1449             :   long n, k;
    1450         476 :   GEN P = cgetg(m+2,t_POL);
    1451         476 :   P[1] = 0;
    1452         476 :   for(n = 0; n < m; n++) gel(P,n+2) = gen_0;
    1453        2205 :   for(n = 0;; n++)
    1454             :   {
    1455        2205 :     k = v * (((3*n - 1) * n) >> 1);
    1456        2205 :     if (k >= m) break;
    1457        1869 :     gel(P, 2+k) = odd(n)? gen_m1: gen_1;
    1458        1869 :     k += n*v; /* v * (3*n + 1) * n / 2 */;
    1459        1869 :     if (k >= m) break;
    1460        1729 :     gel(P, 2+k) = odd(n)? gen_m1: gen_1;
    1461        1729 :   }
    1462         476 :   return RgX_to_ser(P, m+2);
    1463             : }
    1464             : 
    1465             : static GEN
    1466         119 : vecan_eta(GEN eta, long L)
    1467             : {
    1468         119 :   GEN P, eq, divN = gel(eta, 1), RM = gel(eta, 2);
    1469         119 :   long i, ld = lg(divN);
    1470         119 :   P = gen_1; eq = gen_0;
    1471         553 :   for (i = 1; i < ld; ++i)
    1472             :   {
    1473         434 :     GEN m, rm = gel(RM, i);
    1474         434 :     if (!signe(rm)) continue;
    1475         420 :     m = gel(divN, i); eq = addii(eq, mulii(m, rm));
    1476         420 :     P = gmul(P, gpowgs(eta_inflate_ZXn(L, itos(m)), itos(rm)));
    1477             :   }
    1478         119 :   if (!equalis(eq, 24)) pari_err_IMPL("valuation != 1 in lfunetaquo");
    1479         119 :   return gtovec0(P, L);
    1480             : }
    1481             : 
    1482             : /* Check if correct eta quotient. Set canonical form columns */
    1483             : static void
    1484          35 : etaquocheck(GEN eta, GEN *pdivN, GEN *pRM, GEN *pN)
    1485             : {
    1486             :   GEN M, E, divN, RM, N;
    1487             :   long lM, j;
    1488          35 :   if (typ(eta) != t_MAT || lg(eta) != 3 || !RgM_is_ZM(eta))
    1489           0 :     pari_err_TYPE("lfunetaquo", eta);
    1490          35 :   eta = famat_reduce(eta);
    1491          35 :   M = gel(eta, 1); lM = lg(M);
    1492          35 :   E = gel(eta, 2);
    1493          35 :   *pN = N = glcm0(M, NULL);
    1494          35 :   *pdivN = divN = divisors(N);
    1495          35 :   settyp(divN, t_COL);
    1496          35 :   *pRM = RM = zerocol(lg(divN)-1);
    1497         119 :   for (j = 1; j < lM; j++)
    1498             :   {
    1499          84 :     GEN m = gel(M,j), e = gel(E,j);
    1500          84 :     long i = ZV_search(divN, m);
    1501          84 :     gel(RM,i) = e;
    1502             :   }
    1503          35 : }
    1504             : 
    1505             : /* Return etaquotient type:
    1506             :  * -1: nonholomorphic or not on gamma_0(N)
    1507             :  * 0: holomorphic noncuspidal
    1508             :  * 1: cuspidal
    1509             :  * 2: selfdual noncuspidal
    1510             :  * 3: selfdual cuspidal
    1511             :  * Sets conductor, modular weight, canonical matrix */
    1512             : static long
    1513          35 : etaquotype(GEN eta, GEN *pN, long *pw, GEN *pcan)
    1514             : {
    1515             :   GEN divN, RM, S, T, U, N, M;
    1516             :   long ld, i, j, fl;
    1517             : 
    1518          35 :   etaquocheck(eta, &divN, &RM, &N);
    1519          35 :   *pcan = mkmat2(divN, RM);
    1520          35 :   *pw = 0;
    1521          35 :   *pN = gen_1;
    1522             :   /* divN sorted in increasing order, N = last entry, divN[ld-i] = N/divN[i] */
    1523          35 :   ld = lg(divN);
    1524          35 :   S = gen_0; T = gen_0; U = gen_0;
    1525         133 :   for (i = 1; i < ld; ++i)
    1526             :   {
    1527          98 :     GEN m = gel(divN,i), rm = gel(RM,i);
    1528          98 :     if (!signe(rm)) continue;
    1529          84 :     S = addii(S, mulii(m, rm));
    1530          84 :     T = addii(T, rm);
    1531          84 :     U = gadd(U, gdiv(rm, m));
    1532             :   }
    1533          35 :   if (umodiu(S, 24) || umodiu(T, 2)) return -1;
    1534          35 :   *pw = itos(shifti(T,-1));
    1535          35 :   *pN = M = lcmii(N, denom(gdivgs(U, 24)));
    1536         133 :   for (i = 1, fl = 1; i < ld; i++)
    1537             :   {
    1538          98 :     GEN m = gel(divN, i), s = mulii(gel(RM,i), mulii(m,N));
    1539             :     long t;
    1540         448 :     for (j = 1; j < ld; ++j)
    1541         350 :       if (j != i && signe(gel(RM,j)))
    1542             :       {
    1543         210 :         GEN mj = gel(divN, j), nj = gel(divN, ld-j); /* nj = N/mj */
    1544         210 :         s = addii(s, mulii(mulii(gel(RM,j), sqri(gcdii(mj, m))), nj));
    1545             :       }
    1546          98 :     t = signe(s);
    1547          98 :     if (t < 0) return -1;
    1548          98 :     if (t == 0) fl = 0;
    1549             :   }
    1550         133 :   for (i = 1; i < ld; ++i)
    1551             :   {
    1552          98 :     GEN m = gel(divN, i), rm = gel(RM, i);
    1553          98 :     if (!signe(rm)) continue;
    1554          84 :     j = ZV_search(divN, divii(M,m));
    1555          84 :     if (!j || !equalii(rm, gel(RM,j))) return fl;
    1556             :   }
    1557          35 :   return fl+2;
    1558             : }
    1559             : 
    1560             : GEN
    1561          35 : lfunetaquo(GEN eta)
    1562             : {
    1563          35 :   pari_sp ltop = avma;
    1564             :   GEN Ldata, N, can;
    1565             :   long k;
    1566          35 :   switch(etaquotype(eta, &N, &k, &can))
    1567             :   {
    1568          35 :     case 3: break;
    1569           0 :     case 2: pari_err_IMPL("noncuspidal eta quotient");
    1570           0 :     default: pari_err_TYPE("lfunetaquo [non holomorphic]", eta);
    1571             :   }
    1572          35 :   Ldata = mkvecn(6, tag(can, t_LFUN_ETA),
    1573             :                     gen_0, mkvec2(gen_0, gen_1), stoi(k), N, gen_1);
    1574          35 :   return gerepilecopy(ltop, Ldata);
    1575             : }
    1576             : 
    1577             : static GEN
    1578         259 : vecan_qf(GEN Q, long L)
    1579         259 : { return gmul2n(gtovec(qfrep0(Q, utoi(L), 1)), 1); }
    1580             : 
    1581             : long
    1582          98 : qf_iseven(GEN M)
    1583             : {
    1584          98 :   long i, l = lg(M);
    1585         259 :   for (i=1; i<l; i++)
    1586         203 :     if (mpodd(gcoeff(M,i,i))) return 0;
    1587          56 :   return 1;
    1588             : }
    1589             : 
    1590             : GEN
    1591          21 : lfunqf(GEN M, long prec)
    1592             : {
    1593          21 :   pari_sp ltop = avma;
    1594             :   long n, k;
    1595             :   GEN D, d, Mi, cMi, detM, Ldata, poles, res0, res2, eno, dual;
    1596             : 
    1597          21 :   if (typ(M) != t_MAT) pari_err_TYPE("lfunqf", M);
    1598          21 :   if (!RgM_is_ZM(M))   pari_err_TYPE("lfunqf [not integral]", M);
    1599          21 :   n = lg(M)-1;
    1600          21 :   if (odd(n)) pari_err_TYPE("lfunqf [odd dimension]", M);
    1601          21 :   k = n >> 1;
    1602          21 :   M = Q_primpart(M); detM = ZM_det(M);
    1603          21 :   Mi = ZM_inv(M, detM); /* det(M) M^(-1) */
    1604          21 :   if (is_pm1(detM)) cMi = NULL; else Mi = Q_primitive_part(Mi, &cMi);
    1605          21 :   d = cMi? diviiexact(detM, cMi): detM; /* denom(M^(-1)) */
    1606          21 :   if (!qf_iseven(M))
    1607             :   {
    1608          14 :     M = gmul2n(M, 1);
    1609          14 :     d = shifti(d, 1);
    1610          14 :     detM = shifti(detM,n);
    1611             :   }
    1612          21 :   if (!qf_iseven(Mi))
    1613             :   {
    1614          14 :     Mi = gmul2n(Mi, 1);
    1615          14 :     d = shifti(d,1);
    1616             :   }
    1617             :   /* det(Mi) = d^n/det(M), D^2 = det(Mi)/det(M) */
    1618          21 :   D = gdiv(powiu(d,k), detM);
    1619          21 :   if (!issquareall(D, &eno)) eno = gsqrt(D, prec);
    1620          21 :   dual = gequal1(D) ? gen_0: tag(Mi, t_LFUN_QF);
    1621          21 :   res0 = RgX_to_ser(deg1pol_shallow(gen_m2, gen_0, 0), 3);
    1622          21 :   setvalp(res0, -1);
    1623          21 :   res2 = RgX_to_ser(deg1pol_shallow(gmulgs(eno,2), gen_0, 0), 3);
    1624          21 :   setvalp(res2, -1);
    1625          21 :   poles = mkcol2(mkvec2(stoi(k),res2), mkvec2(gen_0,res0));
    1626          21 :   Ldata = mkvecn(7, tag(M, t_LFUN_QF), dual,
    1627             :        mkvec2(gen_0, gen_1), stoi(k), d, eno, poles);
    1628          21 :   return gerepilecopy(ltop, Ldata);
    1629             : }
    1630             : 
    1631             : /********************************************************************/
    1632             : /**  Artin L function, based on a GP script by Charlotte Euvrard   **/
    1633             : /********************************************************************/
    1634             : 
    1635             : static GEN
    1636         119 : artin_repfromgens(GEN G, GEN M)
    1637             : {
    1638         119 :   GEN R, V, ord = gal_get_orders(G), grp = gal_get_group(G);
    1639         119 :   long i, j, k, n = lg(ord)-1, m = lg(grp)-1;
    1640             : 
    1641         119 :   if (lg(M)-1 != n) pari_err_DIM("lfunartin");
    1642         119 :   R = cgetg(m+1, t_VEC);
    1643         119 :   gel(R, 1) = matid(lg(gel(M, 1))-1);
    1644         357 :   for (i = 1, k = 1; i <= n; ++i)
    1645             :   {
    1646         238 :     long c = k*(ord[i] - 1);
    1647         238 :     gel(R, ++k) = gel(M, i);
    1648         238 :     for (j = 2; j <= c; ++j) gel(R, ++k) = gmul(gel(R, j), gel(M, i));
    1649             :   }
    1650         119 :   V = cgetg(m+1, t_VEC);
    1651         119 :   for (i = 1; i <= m; i++) gel(V, gel(grp, i)[1]) = gel(R, i);
    1652         119 :   return V;
    1653             : }
    1654             : 
    1655             : static GEN
    1656          21 : galois_get_conj(GEN G)
    1657             : {
    1658          21 :   GEN grp = gal_get_group(G);
    1659          21 :   long i, k, r = lg(grp)-1;
    1660          21 :   GEN b = zero_F2v(r);
    1661          91 :   for (k = 2; k <= r; ++k)
    1662             :   {
    1663          91 :     GEN g = gel(grp,k);
    1664          91 :     if (!F2v_coeff(b,g[1]) && g[g[1]]==1)
    1665             :     {
    1666          21 :       pari_sp av = avma;
    1667          21 :       GEN F = galoisfixedfield(G, g, 1, -1);
    1668          21 :       if (ZX_sturmpart(F, NULL) > 0) { avma = av; return g; }
    1669           0 :       for (i = 1; i<=r; i++)
    1670             :       {
    1671           0 :         GEN h = gel(grp, i);
    1672           0 :         long t = h[1];
    1673           0 :         while (h[t]!=1) t = h[t];
    1674           0 :         F2v_set(b, h[g[t]]);
    1675             :       }
    1676           0 :       avma = av;
    1677             :     }
    1678             :   }
    1679           0 :   pari_err_BUG("galois_get_conj");
    1680           0 :   return NULL;
    1681             : }
    1682             : 
    1683        7308 : static GEN  cyclotoi(GEN v) { return simplify_shallow(lift_shallow(v)); }
    1684        2681 : static long cyclotos(GEN v) { return gtos(cyclotoi(v)); }
    1685        2660 : static long char_dim(GEN ch) { return cyclotos(gel(ch,1)); }
    1686             : 
    1687             : static GEN
    1688        1330 : artin_gamma(GEN N, GEN G, GEN ch)
    1689             : {
    1690        1330 :   long a, t, d = char_dim(ch);
    1691        1330 :   if (nf_get_r2(N) == 0) return vec01(d, 0);
    1692          21 :   a = galois_get_conj(G)[1];
    1693          21 :   t = cyclotos(gel(ch,a));
    1694          21 :   return vec01((d+t) / 2, (d-t) / 2);
    1695             : }
    1696             : 
    1697             : static long
    1698        3101 : artin_dim(GEN ind, GEN ch)
    1699             : {
    1700        3101 :   long n = lg(ch)-1;
    1701        3101 :   GEN elts = group_elts(ind, n);
    1702        3101 :   long i, d = lg(elts)-1;
    1703        3101 :   GEN s = gen_0;
    1704       17619 :   for(i=1; i<=d; i++)
    1705       14518 :     s = gadd(s, gel(ch, gel(elts,i)[1]));
    1706        3101 :   return gtos(gdivgs(cyclotoi(s), d));
    1707             : }
    1708             : 
    1709             : static GEN
    1710         602 : artin_ind(GEN elts, GEN ch, GEN p)
    1711             : {
    1712         602 :   long i, d = lg(elts)-1;
    1713         602 :   GEN s = gen_0;
    1714        2072 :   for(i=1; i<=d; i++)
    1715        1470 :     s = gadd(s, gel(ch, gmul(gel(elts,i),p)[1]));
    1716         602 :   return gdivgs(s, d);
    1717             : }
    1718             : 
    1719             : static GEN
    1720        2205 : artin_ram(GEN nf, GEN gal, GEN aut, GEN pr, GEN ramg, GEN ch, long d)
    1721             : {
    1722        2205 :   pari_sp av = avma;
    1723             :   long i, v, n;
    1724             :   GEN p, q, V, elts;
    1725        2205 :   if (d==0) return pol_1(0);
    1726         595 :   n = degpol(gal_get_pol(gal));
    1727         595 :   q = p = idealramfrobenius_aut(nf, gal, pr, ramg, aut);
    1728         595 :   elts = group_elts(gel(ramg,2), n);
    1729         595 :   v = fetch_var_higher();
    1730         595 :   V = cgetg(d+3, t_POL);
    1731         595 :   V[1] = evalsigne(1)|evalvarn(v);
    1732         595 :   gel(V,2) = gen_0;
    1733        1197 :   for(i=1; i<=d; i++)
    1734             :   {
    1735         602 :     gel(V,i+2) = gdivgs(artin_ind(elts, ch, q), -i);
    1736         602 :     q = gmul(q, p);
    1737             :   }
    1738         595 :   delete_var();
    1739         595 :   V = RgXn_exp(V,d+1);
    1740         595 :   setvarn(V,0); return gerepileupto(av, ginv(V));
    1741             : }
    1742             : 
    1743             : /* [Artin conductor, vec of [p, Lp]] */
    1744             : static GEN
    1745        1330 : artin_badprimes(GEN N, GEN G, GEN aut, GEN ch)
    1746             : {
    1747        1330 :   pari_sp av = avma;
    1748        1330 :   long i, d = char_dim(ch);
    1749        1330 :   GEN P = gel(absZ_factor(nf_get_disc(N)), 1);
    1750        1330 :   long lP = lg(P);
    1751        1330 :   GEN B = cgetg(lP, t_VEC), C = cgetg(lP, t_VEC);
    1752             : 
    1753        3535 :   for (i = 1; i < lP; ++i)
    1754             :   {
    1755        2205 :     GEN p = gel(P, i), pr = idealprimedec_galois(N, p);
    1756        2205 :     GEN J = idealramgroups_aut(N, G, pr, aut);
    1757        2205 :     GEN G0 = gel(J,2); /* inertia group */
    1758        2205 :     long lJ = lg(J);
    1759        2205 :     long sdec = artin_dim(G0, ch);
    1760        2205 :     long ndec = group_order(G0);
    1761        2205 :     long j, v = ndec * (d - sdec);
    1762        3101 :     for (j = 3; j < lJ; ++j)
    1763             :     {
    1764         896 :       GEN Jj = gel(J, j);
    1765         896 :       long s = artin_dim(Jj, ch);
    1766         896 :       v += group_order(Jj) * (d - s);
    1767             :     }
    1768        2205 :     gel(C, i) = powiu(p, v/ndec);
    1769        2205 :     gel(B, i) = mkvec2(p, artin_ram(N, G, aut, pr, J, ch, sdec));
    1770             :   }
    1771        1330 :   return gerepilecopy(av, mkvec2(ZV_prod(C), B));
    1772             : }
    1773             : 
    1774             : struct dir_artin
    1775             : {
    1776             :   GEN nf, G, V, aut;
    1777             : };
    1778             : 
    1779             : /* p does not divide nf.index */
    1780             : static GEN
    1781       49924 : idealfrobenius_easy(GEN nf, GEN gal, GEN aut, GEN T, GEN p)
    1782             : {
    1783       49924 :   long i, l = lg(aut), f = degpol(T);
    1784       49924 :   GEN D, Dzk, DzkT, DXp, grp = gal_get_group(gal);
    1785       49924 :   pari_sp av = avma;
    1786       49924 :   if (f==1) return gel(grp,1);
    1787       48356 :   Dzk = nf_get_zkprimpart(nf);
    1788       48356 :   D = modii(nf_get_zkden(nf), p);
    1789       48356 :   DzkT = RgV_to_RgM(FqV_red(Dzk, T, p), f);
    1790       48356 :   DXp = RgX_to_RgC(FpX_Frobenius(T, p), f);
    1791       48356 :   if (!equali1(D)) DXp = FpC_Fp_mul(DXp, D, p);
    1792      320488 :   for(i=1; i < l; i++)
    1793             :   {
    1794      320488 :     GEN g = gel(grp,i);
    1795      320488 :     if (perm_order(g)==f)
    1796             :     {
    1797      164010 :       GEN A = FpM_FpC_mul(DzkT, gel(aut,g[1]), p);
    1798      164010 :       if (ZV_equal(A, DXp)) {avma = av; return g; }
    1799             :     }
    1800             :   }
    1801             :   return NULL; /* LCOV_EXCL_LINE */
    1802             : }
    1803             : /* true nf; p divides nf.index, pr/p unramified */
    1804             : static GEN
    1805        1470 : idealfrobenius_hard(GEN nf, GEN gal, GEN aut, GEN pr)
    1806             : {
    1807        1470 :   long i, l = lg(aut), f = pr_get_f(pr);
    1808        1470 :   GEN modpr, p, T, X, Xp, pi, grp = gal_get_group(gal);
    1809        1470 :   pari_sp av = avma;
    1810        1470 :   if (f==1) return gel(grp,1);
    1811        1218 :   pi = pr_get_gen(pr);
    1812        1218 :   modpr = zkmodprinit(nf, pr);
    1813        1218 :   p = modpr_get_p(modpr);
    1814        1218 :   T = modpr_get_T(modpr);
    1815        1218 :   X = modpr_genFq(modpr);
    1816        1218 :   Xp = FpX_Frobenius(T, p);
    1817        8106 :   for (i = 1; i < l; i++)
    1818             :   {
    1819        8106 :     GEN g = gel(grp,i);
    1820        8106 :     if (perm_order(g)==f)
    1821             :     {
    1822        3752 :       GEN S = gel(aut,g[1]);
    1823        3752 :       GEN A = nf_to_Fq(nf, zk_galoisapplymod(nf,X,S,p), modpr);
    1824             :       /* sigma(X) = X^p (mod pr) and sigma(pi) in pr */
    1825        5082 :       if (ZX_equal(A, Xp) && (f == nf_get_degree(nf) ||
    1826        2548 :           ZC_prdvd(zk_galoisapplymod(nf,pi,S,p),pr))) { avma=av; return g; }
    1827             :     }
    1828             :   }
    1829             :   return NULL; /* LCOV_EXCL_LINE */
    1830             : }
    1831             : 
    1832             : static GEN
    1833       51394 : dirartin(void *E, GEN p, long n)
    1834             : {
    1835       51394 :   pari_sp av = avma;
    1836       51394 :   struct dir_artin *d = (struct dir_artin *) E;
    1837       51394 :   GEN nf = d->nf, pr, frob;
    1838             :   /* pick one maximal ideal in the conjugacy class above p */
    1839       51394 :   GEN T = nf_get_pol(nf);
    1840       51394 :   if (!dvdii(nf_get_index(nf), p))
    1841             :   { /* simple case */
    1842       49924 :     GEN F = FpX_factor(T, p), P = gmael(F,1,1);
    1843       49924 :     frob = idealfrobenius_easy(nf, d->G, d->aut, P, p);
    1844             :   }
    1845             :   else
    1846             :   {
    1847        1470 :     pr = idealprimedec_galois(nf,p);
    1848        1470 :     frob = idealfrobenius_hard(nf, d->G, d->aut, pr);
    1849             :   }
    1850       51394 :   avma = av; return RgXn_inv(gel(d->V, frob[1]), n);
    1851             : }
    1852             : 
    1853             : static GEN
    1854        2814 : vecan_artin(GEN an, long L, long prec)
    1855             : {
    1856             :   struct dir_artin d;
    1857        2814 :   GEN A, Sbad = gel(an,5);
    1858        2814 :   long n = itos(gel(an,6)), isreal = lg(an)<8 ? 0: !itos(gel(an,7));
    1859        2814 :   d.nf = gel(an,1); d.G = gel(an,2); d.V = gel(an,3); d.aut = gel(an,4);
    1860        2814 :   A = lift_shallow(direuler_bad(&d, dirartin, gen_2, stoi(L), NULL, Sbad));
    1861        2814 :   A = RgXV_RgV_eval(A, grootsof1(n, prec));
    1862        2814 :   if (isreal) A = real_i(A);
    1863        2814 :   return A;
    1864             : }
    1865             : 
    1866             : static GEN
    1867        2737 : char_expand(GEN conj, GEN ch)
    1868             : {
    1869        2737 :   long i, l = lg(conj);
    1870        2737 :   GEN V = cgetg(l, t_COL);
    1871       30632 :   for (i=1; i<l; i++)
    1872             :   {
    1873       27895 :     long ci = conj[i];
    1874       27895 :     gel(V,i) = gel(ch, ci);
    1875             :   }
    1876        2737 :   return V;
    1877             : }
    1878             : 
    1879             : static GEN
    1880         119 : rep_to_char(GEN R)
    1881             : {
    1882         119 :   long i, l = lg(R);
    1883         119 :   GEN V = cgetg(l, t_COL);
    1884        1281 :   for (i=1; i<l; i++)
    1885        1162 :     gel(V,i) = gtrace(gel(R,i));
    1886         119 :   return V;
    1887             : }
    1888             : 
    1889             : static GEN
    1890        1526 : handle_zeta(long n, GEN ch, long *m)
    1891             : {
    1892             :   GEN c;
    1893        1526 :   long t, i, l = lg(ch);
    1894        1526 :   GEN dim = cyclotoi(vecsum(ch));
    1895        1526 :   if (typ(dim) != t_INT)
    1896           0 :     pari_err_DOMAIN("lfunartin","chi","is not a", strtoGENstr("character"), ch);
    1897        1526 :   t = itos(dim);
    1898        1526 :   if (t < 0 || t % n)
    1899           0 :     pari_err_DOMAIN("lfunartin","chi","is not a", strtoGENstr("character"), ch);
    1900        1526 :   if (t == 0) { *m = 0; return ch; }
    1901         196 :   *m = t / n;
    1902         196 :   c = cgetg(l, t_COL);
    1903        1897 :   for (i=1; i<l; i++)
    1904        1701 :     gel(c,i) = gsubgs(gel(ch,i), *m);
    1905         196 :   return c;
    1906             : }
    1907             : 
    1908             : static int
    1909        6230 : cyclo_is_real(GEN v, GEN ix)
    1910             : {
    1911        6230 :   pari_sp av = avma;
    1912        6230 :   int s = gequal(poleval(lift_shallow(v), ix), v);
    1913        6230 :   avma = av; return s;
    1914             : }
    1915             : 
    1916             : static int
    1917        1330 : char_is_real(GEN ch, GEN mod)
    1918             : {
    1919        1330 :   long i, l = lg(ch);
    1920        1330 :   GEN ix = QXQ_inv(pol_x(varn(mod)), mod);
    1921        6699 :   for (i=1; i<l; i++)
    1922        6230 :     if (!cyclo_is_real(gel(ch,i), ix))
    1923         861 :       return 0;
    1924         469 :   return 1;
    1925             : }
    1926             : 
    1927             : static GEN
    1928        1540 : galois_elts_sorted(GEN gal)
    1929             : {
    1930        1540 :   GEN elts = gal_get_group(gal);
    1931        1540 :   long i, l = lg(elts);
    1932        1540 :   GEN v = cgetg(l, t_VEC);
    1933       16947 :   for(i=1; i<l; i++)
    1934       15407 :     gel(v, mael(elts,i,1)) = gel(elts,i);
    1935        1540 :   return v;
    1936             : }
    1937             : 
    1938             : GEN
    1939        1540 : lfunartin(GEN nf, GEN gal, GEN ch, long o, long bitprec)
    1940             : {
    1941        1540 :   pari_sp av = avma;
    1942        1540 :   GEN bc, V, aut, mod, Ldata = NULL, chx, elts, conj;
    1943             :   long tmult, var, nbc;
    1944        1540 :   nf = checknf(nf);
    1945        1540 :   checkgal(gal);
    1946        1540 :   var = gvar(ch);
    1947        1540 :   if (var == 0) pari_err_PRIORITY("lfunartin",ch,"=",0);
    1948        1540 :   if (var < 0) var = 1;
    1949        1540 :   if (!is_vec_t(typ(ch))) pari_err_TYPE("lfunartin", ch);
    1950        1540 :   elts = galois_elts_sorted(gal);
    1951        1540 :   conj = groupelts_conjclasses(elts, &nbc);
    1952        1540 :   mod = mkpolmod(gen_1, polcyclo(o, var));
    1953        1540 :   if (lg(ch)>1 && typ(gel(ch,1))==t_MAT)
    1954         119 :   {
    1955         119 :     GEN repr = conjclasses_repr(conj, nbc);
    1956         119 :     GEN M = gmul(ch, mod);
    1957         119 :     GEN R = artin_repfromgens(gal, M);
    1958         119 :     chx = rep_to_char(R);
    1959         119 :     ch = shallowextract(chx, repr);
    1960             :   }
    1961             :   else
    1962             :   {
    1963        1421 :     if (nbc != lg(ch)-1) pari_err_DIM("lfunartin");
    1964        1407 :     chx = char_expand(conj, gmul(ch, mod));
    1965             :   }
    1966        1526 :   chx = handle_zeta(nf_get_degree(nf), chx, &tmult);
    1967        1526 :   if (!gequal0(chx))
    1968             :   {
    1969        1330 :     GEN real = char_is_real(chx, gel(mod,1))? gen_0: gen_1;
    1970        1330 :     aut = nfgaloispermtobasis(nf, gal);
    1971        1330 :     V = gmul(char_expand(conj, galoischarpoly(gal, ch, o)), mod);
    1972        1330 :     bc = artin_badprimes(nf, gal, aut, chx);
    1973        3990 :     Ldata = mkvecn(6,
    1974        1330 :       tag(mkcoln(7, nf, gal, V, aut, gel(bc, 2), stoi(o), real), t_LFUN_ARTIN),
    1975        1330 :       real, artin_gamma(nf, gal, chx), gen_1, gel(bc,1), gen_0);
    1976             :   }
    1977        1526 :   if (tmult==0 && Ldata==NULL) pari_err_TYPE("lfunartin",ch);
    1978        1526 :   if (tmult)
    1979             :   {
    1980             :     long i;
    1981         196 :     if (Ldata==NULL) { Ldata = lfunzeta(); tmult--; }
    1982         196 :     for(i=1; i<=tmult; i++)
    1983           0 :       Ldata = lfunmul(Ldata, gen_1, bitprec);
    1984             :   }
    1985        1526 :   return gerepilecopy(av, Ldata);
    1986             : }
    1987             : 
    1988             : /********************************************************************/
    1989             : /**                    High-level Constructors                     **/
    1990             : /********************************************************************/
    1991             : enum { t_LFUNMISC_POL, t_LFUNMISC_CHIQUAD, t_LFUNMISC_CHICONREY,
    1992             :        t_LFUNMISC_CHIGEN, t_LFUNMISC_ELLINIT, t_LFUNMISC_ETAQUO };
    1993             : static long
    1994        2240 : lfundatatype(GEN data)
    1995             : {
    1996             :   long l;
    1997        2240 :   switch(typ(data))
    1998             :   {
    1999         644 :     case t_INT: return t_LFUNMISC_CHIQUAD;
    2000           7 :     case t_INTMOD: return t_LFUNMISC_CHICONREY;
    2001         399 :     case t_POL: return t_LFUNMISC_POL;
    2002             :     case t_VEC:
    2003        1190 :       if (checknf_i(data)) return t_LFUNMISC_POL;
    2004        1078 :       l = lg(data);
    2005        1078 :       if (l == 17) return t_LFUNMISC_ELLINIT;
    2006         301 :       if (l == 3 && typ(gel(data,1)) == t_VEC) return t_LFUNMISC_CHIGEN;
    2007           0 :       break;
    2008             :   }
    2009           0 :   return -1;
    2010             : }
    2011             : static GEN
    2012       21397 : lfunmisc_to_ldata_i(GEN ldata, long shallow)
    2013             : {
    2014             :   long lx;
    2015       21397 :   if (is_linit(ldata)) ldata = linit_get_ldata(ldata);
    2016       21397 :   lx = lg(ldata);
    2017       21397 :   if (typ(ldata)==t_VEC && (lx == 7 || lx == 8) && is_tagged(ldata))
    2018             :   {
    2019       19157 :     if (!shallow) ldata = gcopy(ldata);
    2020       19157 :     checkldata(ldata); return ldata;
    2021             :   }
    2022        2240 :   switch (lfundatatype(ldata))
    2023             :   {
    2024         511 :     case t_LFUNMISC_POL: return lfunzetak(ldata);
    2025         644 :     case t_LFUNMISC_CHIQUAD: return lfunchiquad(ldata);
    2026             :     case t_LFUNMISC_CHICONREY:
    2027             :     {
    2028           7 :       GEN G = znstar0(gel(ldata,1), 1);
    2029           7 :       return lfunchiZ(G, gel(ldata,2));
    2030             :     }
    2031         301 :     case t_LFUNMISC_CHIGEN: return lfunchigen(gel(ldata,1), gel(ldata,2));
    2032         777 :     case t_LFUNMISC_ELLINIT: return lfunell(ldata);
    2033             :   }
    2034           0 :   pari_err_TYPE("lfunmisc_to_ldata",ldata);
    2035             :   return NULL; /* LCOV_EXCL_LINE */
    2036             : }
    2037             : 
    2038             : GEN
    2039         217 : lfunmisc_to_ldata(GEN ldata)
    2040         217 : { return lfunmisc_to_ldata_i(ldata, 0); }
    2041             : 
    2042             : GEN
    2043       21180 : lfunmisc_to_ldata_shallow(GEN ldata)
    2044       21180 : { return lfunmisc_to_ldata_i(ldata, 1); }
    2045             : 
    2046             : /********************************************************************/
    2047             : /**                    High-level an expansion                     **/
    2048             : /********************************************************************/
    2049             : /* van is the output of ldata_get_an: return a_1,...a_L at precision prec */
    2050             : GEN
    2051       10409 : ldata_vecan(GEN van, long L, long prec)
    2052             : {
    2053       10409 :   GEN an = gel(van, 2);
    2054       10409 :   long t = mael(van,1,1);
    2055             :   pari_timer ti;
    2056       10409 :   if (DEBUGLEVEL >= 1)
    2057           0 :     err_printf("Lfun: computing %ld coeffs, prec %ld, type %ld\n", L, prec, t);
    2058       10409 :   if (DEBUGLEVEL >= 2) timer_start(&ti);
    2059       10409 :   switch (t)
    2060             :   {
    2061             :     long n;
    2062             :     case t_LFUN_GENERIC:
    2063         665 :       an = vecan_closure(an, L, prec);
    2064         658 :       n = lg(an)-1;
    2065         658 :       if (n < L)
    2066           7 :         pari_warn(warner, "#an = %ld < %ld, results may be imprecise", n, L);
    2067         658 :       break;
    2068        1673 :     case t_LFUN_ZETA: an = const_vecsmall(L, 1); break;
    2069        1232 :     case t_LFUN_NF:  an = dirzetak(an, stoi(L)); break;
    2070         826 :     case t_LFUN_ELL: an = ellan(an, L); break;
    2071         574 :     case t_LFUN_KRONECKER: an = vecan_Kronecker(an, L); break;
    2072         294 :     case t_LFUN_CHIZ: an = vecan_chiZ(an, L, prec); break;
    2073         735 :     case t_LFUN_CHIGEN: an = vecan_chigen(an, L, prec); break;
    2074        2814 :     case t_LFUN_ARTIN: an = vecan_artin(an, L, prec); break;
    2075         119 :     case t_LFUN_ETA: an = vecan_eta(an, L); break;
    2076         259 :     case t_LFUN_QF: an = vecan_qf(an, L); break;
    2077         609 :     case t_LFUN_DIV: an = vecan_div(an, L, prec); break;
    2078         161 :     case t_LFUN_MUL: an = vecan_mul(an, L, prec); break;
    2079           0 :     case t_LFUN_CONJ: an = vecan_conj(an, L, prec); break;
    2080          56 :     case t_LFUN_SYMSQ_ELL: an = vecan_ellsymsq(an, L); break;
    2081          28 :     case t_LFUN_GENUS2: an = vecan_genus2(an, L); break;
    2082         364 :     case t_LFUN_MFCLOS: an = mfcoefs(an,L,1) + 1; /* skip a_0 */
    2083         364 :                         an[0] = evaltyp(t_VEC)|evallg(L+1); break;
    2084           0 :     default: pari_err_TYPE("ldata_vecan", van);
    2085             :   }
    2086       10402 :   if (DEBUGLEVEL >= 2) timer_printf(&ti, "ldata_vecan");
    2087       10402 :   return an;
    2088             : }

Generated by: LCOV version 1.11