Code coverage tests

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

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

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

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

Generated by: LCOV version 1.13