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 - hgm.c (source / functions) Hit Total Coverage
Test: PARI/GP v2.18.0 lcov report (development 29792-e7296d7d08) Lines: 1219 1258 96.9 %
Date: 2024-12-05 09:09:40 Functions: 139 139 100.0 %
Legend: Lines: hit not hit

          Line data    Source code
       1             : /* Copyright (C) 2018  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             : /*                   HYPERGEOMETRIC MOTIVES                     */
      16             : /****************************************************************/
      17             : #include "pari.h"
      18             : #include "paripriv.h"
      19             : 
      20             : #define DEBUGLEVEL DEBUGLEVEL_hgm
      21             : static THREAD GEN HGMDATA2, HGMDATA3;
      22             : 
      23             : void
      24      326281 : pari_init_hgm(void) { HGMDATA2 = HGMDATA3 = NULL; }
      25             : 
      26             : void
      27      321516 : pari_close_hgm(void)
      28      321516 : { guncloneNULL_deep(HGMDATA2); guncloneNULL_deep(HGMDATA3); }
      29             : 
      30             : /* For now: only HGM defined over Q. Can be given as
      31             :  * [alpha] two lists (alpha_i), (beta_j) of the same length of
      32             :  *      rational numbers in [0,1[ with the cyclotomic property.
      33             :  * [cyclo] two lists (d_i), (e_j) of positive integers with
      34             :  *      sum_i phi(d_i) = sum_j phi(e_j).
      35             :  * [gamma] a list of (g_i) such that P=prod_i(x^i-1)^{g_i}
      36             :  *
      37             :  * Initialization computes all HGM data independent of t (hgm) */
      38             : /***************************************************************/
      39             : /*            PART I: functions independent of t, p            */
      40             : /***************************************************************/
      41             : 
      42             : static GEN
      43        1064 : RgV_addhalf(GEN x) { pari_APPLY_same(gadd(ghalf, gel(x,i))); }
      44             : static GEN
      45      979433 : RgV_mirror(GEN x) { pari_APPLY_same(gsubsg(1, gel(x,i))); }
      46             : 
      47             : static GEN
      48      104146 : hodge(GEN val, GEN vbe, long *pTT)
      49             : {
      50      104146 :   long r = lg(val) - 1, T, s, k;
      51      104146 :   GEN H, w = indexsort( shallowconcat(val, RgV_mirror(vbe)) );
      52     1854720 :   for (k = 1, s = T = 0; k <= 2*r; k++)
      53     1750574 :     if (w[k] > r) { if (--s < T) T = s; } else s++;
      54      104146 :   H = zero_zv(r + 1 - T);
      55     1854720 :   for (k = 1, s = 0; k <= 2*r; k++) /* (x^-T * Hodge) += x^s, s - T >= 0 */
      56     1750574 :     if (w[k] > r) s--; else { H[s + 1 - T]++; s++; }
      57      104146 :   *pTT = -T; return zv_to_zx(H, evalvarn(0));
      58             : }
      59             : 
      60             : /* Conversion from [val] to cyclotomic factorization: t_VECSMALL E encodes
      61             :  * prod_i polcyclo(i)^E[i] */
      62             : static GEN
      63        1414 : al2cyE(GEN v)
      64             : {
      65             :   GEN dv, E, F, w;
      66        1414 :   long N, i, j, l = lg(v);
      67        1414 :   if (l == 1) return cgetg(1, t_VECSMALL);
      68        1246 :   w = Q_remove_denom(v, &dv);
      69        1246 :   if (!dv) return mkvecsmall(l-1);
      70        1001 :   N = itou(dv); w = ZV_to_Flv(w, N); vecsmall_sort(w);
      71        1001 :   E = zero_zv(N); /* to record polcyclo(i)^E[i] */
      72        1001 :   F = cgetg(l, t_VECSMALL); /* to check input */
      73        4949 :   for (i = j = 1; i < l; i++)
      74             :   {
      75        3948 :     long k, m = w[i];
      76        3948 :     if (!m) { E[1]++; F[j++] = 0; }
      77        3563 :     else if (N % m == 0)
      78             :     {
      79        1904 :       long d = N/m, km;
      80        1904 :       E[d]++;
      81        9037 :       for (k = 1, km = m; k <= d; k++, km += m)
      82        7133 :         if (ugcd(k,d) == 1) F[j++] = km;
      83             :     }
      84             :   }
      85        1001 :   setlg(F,j); vecsmall_sort(F);
      86        1001 :   return gequal(w,F)? E: NULL;
      87             : }
      88             : 
      89             : static int
      90         707 : cyE_intersect(GEN a, GEN b)
      91             : {
      92         707 :   long i, l = minss(lg(a), lg(b));
      93        2695 :   for (i = 1; i < l; i++) if (a[i] && b[i]) return 1;
      94         700 :   return 0;
      95             : }
      96             : static GEN
      97         707 : get_CYCLOE(GEN val, GEN vbe)
      98             : {
      99         707 :   GEN Ea = al2cyE(val), Eb = al2cyE(vbe);
     100         707 :   if (!Ea || !Eb || cyE_intersect(Ea, Eb))
     101           7 :     pari_err_TYPE("hgminit [not a Q motive]", mkvec2(val,vbe));
     102         700 :   return mkvec2(Ea, Eb);
     103             : }
     104             : 
     105             : /* R[n/d] += mu(d) * r, for all d | n */
     106             : static void
     107       10549 : moebiusadd(GEN R, long n, long r)
     108             : {
     109       10549 :   if (r) {
     110        1785 :     GEN D = divisorsu_moebius(gel(factoru(n), 1));
     111        1785 :     long j, l = lg(D);
     112        1785 :     R[n] += r; /* d = 1 */
     113        3696 :     for (j = 2; j < l; j++) { long d = D[j]; R[n / labs(d)] += d < 0? -r: r; }
     114             :   }
     115       10549 : }
     116             : /* Use Phi_n(x) = Prod_{d|n} (x^d-1)^mu(n/d); Ea/Eb in terms of (x^i-1) */
     117             : static GEN
     118         700 : get_VPOLGA(GEN E)
     119             : {
     120         700 :   GEN Ea = gel(E,1), Eb = gel(E,2), R;
     121         700 :   long i, la = lg(Ea), lb = lg(Eb), n = maxuu(la, lb) - 1;
     122             :   pari_sp av;
     123         700 :   R = zero_zv(n); av = avma;
     124        7770 :   for (i = 1; i < la; i++) moebiusadd(R, i, Ea[i]);
     125        4179 :   for (i = 1; i < lb; i++) moebiusadd(R, i, -Eb[i]);
     126        5733 :   for (i = n; i; i--) if (R[i]) break;
     127         700 :   setlg(R,i+1); set_avma(av); return R;
     128             : }
     129             : 
     130             : /* disc(polcyclo(n)) modulo squares */
     131             : static GEN
     132         126 : cyclodiscmodsq(ulong n)
     133             : {
     134             :   long e, s;
     135             :   ulong p;
     136             : 
     137         126 :   if (n <= 2) return gen_1;
     138         126 :   if ((n & 3) == 2) n >>= 1;
     139         126 :   e = uisprimepower(n, &p);
     140         126 :   if (!e) return gen_1;
     141         112 :   if (p == 2) return e == 2? gen_m1: gen_1;
     142          63 :   s = odd(e)? p: 1;
     143          63 :   if ((p & 3) == 3) s = -s;
     144          63 :   return stoi(s);
     145             : }
     146             : /* E encodes prod polcyclo(i)^E(i) */
     147             : static GEN
     148         231 : discprod(GEN E)
     149             : {
     150         231 :   long i, j, l = lg(E);
     151         231 :   GEN P = cgetg(l, t_VEC);
     152        1148 :   for (i = 3, j = 1; i < l; i++) /* can skip polcyclo(1 or 2) */
     153         917 :     if (odd(E[i])) gel(P,j++) = cyclodiscmodsq(i);
     154         231 :   setlg(P,j); return ZV_prod(P);
     155             : }
     156             : 
     157             : static GEN
     158        1218 : QV_normalize(GEN a)
     159             : {
     160        1218 :   long l = lg(a), i;
     161        1218 :   GEN v = cgetg(l, t_VEC);
     162        5978 :   for (i = 1; i < l; i++)
     163             :   {
     164        4774 :     GEN c = gel(a, i);
     165        4774 :     if (!is_rational_t(typ(c)))
     166          14 :       pari_err_TYPE("hgminit [not rational params]",c);
     167        4760 :     gel(v, i) = gfrac(c);
     168             :   }
     169        1204 :   return sort(v);
     170             : }
     171             : 
     172             : static GEN
     173        1365 : mangoldtexp(long n)
     174             : {
     175        1365 :   GEN D = divisorsu_moebius(gel(factoru(n),1)), P = gen_1;
     176        1365 :   long i, l = lg(D);
     177        4151 :   for (i = 1; i < l; i++)
     178             :   { /* prod (n/d) ^ (n/d mu(d)) */
     179        2786 :     long mud = D[i], nd = n / mud; /* mud = mu(d) * d */
     180        2786 :     P = gmul(P, powis(utoi(labs(nd)), nd));
     181             :   }
     182        1365 :   return P;
     183             : }
     184             : static GEN
     185         938 : E2exp(GEN E)
     186             : {
     187         938 :   long n, r, l = lg(E);
     188         938 :   GEN P = gen_1;
     189        9877 :   for (n = 1; n < l; n++)
     190        8939 :     if ((r = E[n])) P = gmul(P, gpowgs(mangoldtexp(n), r));
     191         938 :   return P;
     192             : }
     193             : static long
     194         973 : get_b1(GEN CYCLOE) { return gel(CYCLOE,2)[1]; }
     195             : 
     196             : static GEN get_u(GEN al, GEN be, GEN CYCLOE, GEN VPOLGA, long DEG, long WT);
     197             : /* [a,b] = alpha, beta */
     198             : static GEN
     199         511 : initab(GEN a, GEN b)
     200             : {
     201             :   GEN VALNUM, VALDEN, VBENUM, VBEDEN, CYCLOE, SIGNPAR;
     202             :   GEN BAD, HODGE, VPOLGA, res;
     203         511 :   long j, WT, TT, OFFMPOL, l = lg(a), DEG = l-1, SWAP = 0;
     204             : 
     205         511 :   if (lg(b) != l) pari_err_TYPE("hgminit [#al != #be]", mkvec2(a,b));
     206         511 :   if (l == 1) pari_err_TYPE("hgminit [#al = 0]", mkvec2(a,b));
     207         497 :   a = QV_normalize(a);
     208         490 :   b = QV_normalize(b);
     209         483 :   if (isintzero(gel(a,1)))
     210             :   {
     211          49 :     if (isintzero(gel(b,1)))
     212           7 :       pari_err_TYPE("hgminit [nonempty intersection]", mkvec2(a,b));
     213          42 :     swap(a, b); SWAP = 1;
     214             :   }
     215         476 :   VALNUM = cgetg(l, t_VECSMALL);
     216         476 :   VALDEN = cgetg(l, t_VECSMALL);
     217         476 :   VBENUM = cgetg(l, t_VECSMALL);
     218         476 :   VBEDEN = cgetg(l, t_VECSMALL);
     219        2478 :   for (j = 1; j < l; j++)
     220             :   {
     221        2002 :     Qtoss(gel(a,j), &VALNUM[j], &VALDEN[j]);
     222        2002 :     Qtoss(gel(b,j), &VBENUM[j], &VBEDEN[j]);
     223             :   }
     224         476 :   BAD = negi(lcmii(ZV_lcm(zv_to_ZV(VALDEN)), ZV_lcm(zv_to_ZV(VBEDEN))));
     225         476 :   HODGE = hodge(a, b, &TT);
     226         476 :   WT = degpol(HODGE);
     227         476 :   CYCLOE = get_CYCLOE(a, b);
     228         469 :   VPOLGA = get_VPOLGA(CYCLOE);
     229         469 :   SIGNPAR = odd(WT)? gen_1: discprod(gel(CYCLOE, odd(DEG)? 2: 1));
     230         469 :   OFFMPOL = (get_b1(CYCLOE) - zv_sum(VPOLGA)) >> 1;
     231         469 :   res = cgetg(13, t_VEC);
     232         469 :   gel(res, 1) = a;
     233         469 :   gel(res, 2) = b;
     234         469 :   gel(res, 3) = CYCLOE;
     235         469 :   gel(res, 4) = VBENUM;
     236         469 :   gel(res, 5) = VBEDEN;
     237         469 :   gel(res, 6) = gdiv(E2exp(gel(CYCLOE,1)), E2exp(gel(CYCLOE,2))); /*MVALUE*/
     238         469 :   gel(res, 7) = VPOLGA;
     239         469 :   gel(res, 8) = BAD;
     240         469 :   gel(res, 9) = HODGE;
     241         469 :   gel(res, 10) = get_u(a, b, CYCLOE, VPOLGA, DEG, WT);
     242         469 :   gel(res, 11) = SIGNPAR;
     243         469 :   gel(res, 12) = mkvecsmall3(OFFMPOL, TT, SWAP);
     244         469 :   return res;
     245             : }
     246             : static int
     247        1295 : checkhgm(GEN v)
     248             : {
     249        1274 :   return typ(v) == t_VEC && lg(v) == 13
     250        2569 :       && typ(gel(v,12)) == t_VECSMALL && lg(gel(v,12)) == 4;
     251             : }
     252             : 
     253             : /* accessors */
     254      141826 : static GEN hgm_get_VAL(GEN v) { return gel(v, 1); }
     255          35 : static GEN hgm_get_VBE(GEN v) { return gel(v, 2); }
     256       12446 : static GEN hgm_get_CYCLOE(GEN v) { return gel(v, 3); }
     257       37597 : static GEN hgm_get_VBENUM(GEN v) { return gel(v, 4); }
     258       37597 : static GEN hgm_get_VBEDEN(GEN v) { return gel(v, 5); }
     259         693 : static GEN hgm_get_MVALUE(GEN v) { return gel(v, 6); }
     260      698752 : static GEN hgm_get_VPOLGA(GEN v) { return gel(v, 7); }
     261       37846 : static GEN hgm_get_BAD(GEN v) { return gel(v, 8); }
     262      109812 : static GEN hgm_get_HODGE(GEN v) { return gel(v, 9); }
     263        2940 : static GEN hgm_get_U(GEN v) { return gmael(v, 10,1); }
     264          63 : static GEN hgm_get_U0(GEN v) { return gmael(v, 10,2); }
     265       28630 : static GEN hgm_get_SIGNPAR(GEN v) { return gel(v, 11); }
     266      141793 : static long hgm_get_DEG(GEN v) { return lg(hgm_get_VAL(v))-1; }
     267      652499 : static long hgm_get_OFFMPOL(GEN v) { return gel(v, 12)[1]; }
     268       87838 : static long hgm_get_TT(GEN v) { return gel(v, 12)[2]; }
     269         686 : static long hgm_get_SWAP(GEN v) { return gel(v, 12)[3]; }
     270      109449 : static long hgm_get_WT(GEN v) { return degpol(hgm_get_HODGE(v)); }
     271             : 
     272             : /***************************************************************/
     273             : /*                 PART II: p-adic functions                   */
     274             : /***************************************************************/
     275             : 
     276             : /* Compute p^k u_{pk+r} for 0 <= r < p and 0 <= k < ell */
     277             : static GEN
     278        2155 : compu(long ell, long p, long D)
     279             : {
     280        2155 :   GEN v = cgetg(p + 1, t_VEC), w = cgetg(ell + 2, t_VEC);
     281        2155 :   long s, t, k, i, L = ell * p;
     282             : 
     283        2155 :   gel(v,1) = gel(v,2) = cvtop(gen_1, stoi(p), D);
     284       59342 :   for (i = 2; i < p; i++) gel(v, i+1) = gdivgu(gel(v, i), i);
     285        2155 :   gel(w, 1) = shallowcopy(v);
     286      261235 :   for (k = s = 1, t = p; k <= L; k++)
     287             :   {
     288      259082 :     gel(v, s) = gdivgu(gadd(gel(v, t), gel(v, s)), k + p - 1);
     289      259080 :     if (++s > p) { s = 1; gel(w, k/p + 1) = shallowcopy(v); }
     290      259080 :     if (++t > p) t = 1;
     291             :   }
     292        2153 :   return w;
     293             : }
     294             : 
     295             : /* [binom(x + k - 1, k) k! p^k, k = 1..N], x a t_PADIC */
     296             : static GEN
     297      116953 : binomfact(GEN x, long N)
     298             : {
     299      116953 :   GEN z, v = cgetg(N+1, t_VEC);
     300             :   long i;
     301      116953 :   gel(v,1) = z = leafcopy(x);
     302      116952 :   setvalp(z, valp(z)+1);
     303      625628 :   for (i = 1; i < N; i++)
     304             :   {
     305      508692 :     gel(v, i+1) = z = gmul(z, gaddgs(x, i));
     306      508641 :     setvalp(z, valp(z) + 1);
     307             :   }
     308      116936 :   return v;
     309             : }
     310             : /* gamma(m /p^f-1) mod p^dfp, m = m0 (mod p), x = (m\p + m0*p^(f-1)) / p^f-1 */
     311             : static ulong
     312      116953 : sumbinom(GEN v, GEN x, long m0, long N, ulong pD)
     313             : {
     314      116953 :   pari_sp av = avma;
     315      116953 :   GEN C = binomfact(x, N), S = gmael(v, 1, m0+1); /* k = 0 */
     316             :   long k;
     317      742533 :   for (k = 1; k <= N; k++) S = gadd(S, gmul(gel(C,k), gmael(v, k+1, m0+1)));
     318      116943 :   return gc_ulong(av, Rg_to_Fl(S, pD));
     319             : }
     320             : /* x * y mod p^D, y a p-unit */
     321             : static GEN
     322      203900 : umultop(ulong x, ulong y, GEN gp, GEN gpD, long D)
     323             : {
     324      203900 :   ulong pD = gpD[2];
     325             :   long v;
     326             :   GEN z;
     327      203900 :   if (!x) return zeropadic_shallow(gp, D);
     328      203900 :   v = u_lvalrem(x, gp[2], &x);
     329      203900 :   if (x >= pD) x %= pD;
     330      203900 :   z = cgetg(5, t_PADIC);
     331      203900 :   z[1] = evalprecp(D) | evalvalp(v);
     332      203899 :   gel(z,2) = gp;
     333      203899 :   gel(z,3) = gpD;
     334      203899 :   gel(z,4) = utoi( Fl_mul(x, y, pD) ); return z;
     335             : }
     336             : 
     337             : /* Compute all gamma(m/(p^f-1)) mod p^D for 0 <= m < p^f-1 */
     338             : static GEN
     339        2162 : precomp(ulong p, long f, long D)
     340             : {
     341             :   pari_sp av, av2;
     342             :   GEN vall, gp, gpD, vres, x, l2;
     343             :   ulong m, m0, m1, M, q, v, vpo, vva, vpowv, pf1;
     344        2162 :   ulong N = ceildivuu(D * p, p-1), pf = upowuu(p,f), pD = upowuu(p, D), pM;
     345        2162 :   ulong ga12, S, iQ, Q = pf - 1;
     346             : 
     347        2162 :   if (!pD || !pf) pari_err_OVERFLOW("hgmprecomp");
     348        2155 :   iQ = D <= f? pD - 1: Fl_inv(Q, pD);
     349        2155 :   vres = cgetg(pf, t_VECSMALL); vres[1] = 1; av = avma;
     350        2155 :   gp = utoipos(p); gpD = utoipos(pD);
     351        2155 :   vall = compu(N, p, D);
     352        2155 :   if (p == 2)
     353             :   {
     354         133 :     if (f == 1) return gc_const(av, vres);
     355          14 :     av2 = avma;
     356          70 :     for (m = 1; m < Q; m++, set_avma(av2))
     357             :     {
     358          56 :       m0 = m&1, m1 = m >> 1;
     359          56 :       x = umultop(m1 + (m0 << (f - 1)), iQ, gp, gpD, D);
     360          56 :       vres[m+1] = sumbinom(vall, x, m0, N, pD);
     361             :     }
     362          14 :     return gc_const(av, vres);
     363             :   }
     364        2022 :   q = Q >> 1; pf1 = pf / p;
     365        2022 :   vva = z_lval(Q, 2); vpowv = 1 << vva;
     366        2022 :   av2 = avma;
     367      104483 :   for (m = 1; m <= q; m += 2, set_avma(av2))
     368             :   {
     369      102461 :     if (f > 1) { m0 = m%p; m1 = m / p; M = m1 + m0 * pf1; } else M = m0 = m;
     370      102461 :     x = umultop(M, iQ, gp, gpD, D);
     371      102460 :     vres[m+1] = S = sumbinom(vall, x, m0, N, pD);
     372      102457 :     if (!odd(m0)) S = pD - S;
     373      102457 :     vres[pf - m] = Fl_inv(S, pD);
     374             :   }
     375       16460 :   for (m = vpowv; m <= q; m += 2*vpowv, set_avma(av2))
     376             :   {
     377       14438 :     if (f > 1) { m0 = m%p; m1 = m / p; M = m1 + m0 * pf1; } else M = m0 = m;
     378       14438 :     x = umultop(M, iQ, gp, gpD, D);
     379       14438 :     vres[m+1] = S = sumbinom(vall, x, m0, N, pD);
     380       14438 :     if (!odd(m0)) S = pD - S;
     381       14438 :     vres[pf - m] = Fl_inv(S, pD);
     382             :   }
     383        2022 :   l2 = gmulgu(Qp_log(cvtop(gen_2, gp, D)), p - 1);
     384        2022 :   ga12 = Fl_inv(Rg_to_Fl(Qp_gamma(cvtop(ghalf, gp, D)), pD), pD);
     385        2022 :   pM = maxuu(pD, pf); av2 = avma;
     386       10637 :   for (v = 1, vpo = 2; vpo < Q; v++, vpo <<= 1)
     387             :   {
     388        8615 :     if (v == vva) continue;
     389       93797 :     for (m = vpo; m <= q; m += 2*vpo, set_avma(av2))
     390             :     {
     391       86947 :       ulong ms2 = (m >> 1) + 1;
     392       86947 :       ulong t = Fl_mul(Fl_mul(vres[ms2], vres[ms2+q], pD), ga12, pD);
     393             :       GEN mp;
     394             :       /* - m  = -m1 * p + m0, 1 <= m0 <= p */
     395       86947 :       m0 = p - m%p; m1 = m/p + 1; M = pM + m1 - m0 * pf1;
     396       86947 :       mp = umultop(M, iQ, gp, gpD, D);
     397       86945 :       S = Rg_to_Fl(gmul(Qp_exp(gmul(mp, l2)), shifti(utoi(t), m0-1)), pD);
     398       86946 :       vres[m+1] = S;
     399       86946 :       if (odd(m0)) S = pD - S;
     400       86946 :       vres[pf - m] = Fl_inv(S, pD);
     401             :     }
     402             :   }
     403        2022 :   return gc_const(av, vres);
     404             : }
     405             : 
     406             : static long
     407        3661 : get_pad(ulong p)
     408             : {
     409        3661 :   switch(p) {
     410         182 :     case 2: return 18;
     411         399 :     case 3: return 11;
     412         322 :     case 5: return  8;
     413             :   }
     414        2758 :   if (p <= 37) return 6;
     415         819 :   if (p <= 251) return 4;
     416           7 :   if (p <= 65521) return 2;
     417           7 :   return 1;
     418             : }
     419             : 
     420             : static GEN
     421         609 : Flv_red(GEN z, ulong p)
     422             : {
     423         609 :   long i, l = lg(z);
     424         609 :   GEN x = cgetg(l, t_VECSMALL);
     425      248780 :   for (i=1; i<l; i++) x[i] = uel(z,i)%p;
     426         609 :   return x;
     427             : }
     428             : 
     429             : static GEN
     430          97 : block0(long l)
     431             : {
     432          97 :   GEN v = cgetg_block(l, t_VEC);
     433             :   long i;
     434       89060 :   for (i = 1; i < l; i++) gel(v,i) = gen_0;
     435          97 :   return v;
     436             : }
     437             : 
     438             : /* f > 0, p prime, need result mod p^dfp */
     439             : static GEN
     440        3668 : doprecomp(ulong p, long f, long dfp)
     441             : {
     442             :   GEN t, T;
     443             :   long F; /* can compute cache mod p^F < 2^BIL */
     444             :   ulong m, n, q;
     445        3668 :   if (f > 3 || (F = get_pad(p)) < dfp) return precomp(p, f, dfp);
     446        3633 :   if (f == 3)
     447             :   {
     448          35 :     T = HGMDATA3; if (!T) T = HGMDATA3 = block0(100);
     449          35 :     if (p >= (ulong)lg(T)) return precomp(p, f, dfp);
     450          35 :     t = gel(T,p);
     451          35 :     if (typ(t) == t_INT) t = gel(T,p) = gclone(precomp(p, f, F));
     452          35 :     return (F == dfp)? t: Flv_red(t, upowuu(p, dfp));
     453             :   }
     454             :   /* f <= 2, dfp <= F */
     455        3598 :   T = HGMDATA2; if (!T) T = HGMDATA2 = block0(1000);
     456        3598 :   if (p >= (ulong)lg(T)) return precomp(p, f, dfp);
     457        3598 :   t = gel(T,p);
     458        3598 :   if (typ(t) == t_INT)
     459             :   {
     460        2099 :     if (f == 1) return precomp(p, f, dfp);
     461         133 :     t = gel(T,p) = gclone(precomp(p, f, F));
     462             :   }
     463             :   /* t = precomp(p, 2, F) = HGMDATA2[p] */
     464        1632 :   if (f == 2)
     465         651 :     return (F == dfp)? t: Flv_red(t, upowuu(p, dfp));
     466             :   /* f = 1 */
     467         981 :   T = cgetg(p, t_VECSMALL); q = upowuu(p, dfp);
     468       10909 :   for (m = n = 1; m < p; m++, n += p+1) T[m] = uel(t, n) % q;
     469         981 :   return T;
     470             : }
     471             : 
     472             : /* W[N] = teich(N ^ (N * sign(VPOLGA[N]))) mod p^2 */
     473             : static GEN
     474       33936 : teichmodp2(GEN VPOLGA, ulong p, ulong p2)
     475             : {
     476             :   ulong pN;
     477             :   long c, l, N;
     478       33936 :   GEN W = cgetg_copy(VPOLGA,&l); /* W[N] = teich(N)^N */
     479      142863 :   for (N = 1, pN = p; N < l; N++, pN += p)
     480      108931 :     if ((c = VPOLGA[N]))
     481             :     {
     482       82251 :       long N0; (void)z_lvalrem(N, p, &N0);
     483       82248 :       if (c < 0) N0 = Fl_inv(N0 % p, p);
     484       82249 :       W[N] = Fl_powu(N0, pN, p2);
     485             :     }
     486       33932 :   return W;
     487             : }
     488             : /* GF[i] = i! mod p^2 */
     489             : static GEN
     490       33928 : get_GF(ulong p, ulong p2)
     491             : {
     492       33928 :   GEN F = cgetg(p, t_VECSMALL);
     493             :   ulong i;
     494    33570292 :   F[1] = 1; F[2] = 2; for (i = 3; i < p; i++) F[i] = Fl_mul(F[i-1], i, p2);
     495       31438 :   return F;
     496             : }
     497             : 
     498             : /* p odd, return v s.t. v[i] = 1/i mod p, i < p/2 */
     499             : static GEN
     500       33934 : Flv_inv_p2(ulong p)
     501             : {
     502       33934 :   long i, g = pgener_Fl(p), ph = (p >> 1);
     503       33934 :   GEN v = cgetg(ph + 1, t_VECSMALL);
     504       33935 :   pari_sp av = avma;
     505       33935 :   GEN w = cgetg(ph, t_VECSMALL);
     506       34893 :   w[1] = g;
     507    16914099 :   for (i = 2; i < ph; i++) w[i] = Fl_mul(w[i-1], g, p); /* w[i]=g^i */
     508    17371941 :   for (i = 1; i < ph; i++)
     509             :   {
     510    17338018 :     long x = w[i], y = w[ph - i]; /* g^(-i) = - g^(ph - i) */
     511    17338018 :     if (x > ph) v[p - x] = y; else v[x] = p - y; /* 1/(-x) = y, 1/x = -y */
     512             :   }
     513       33923 :   v[1] = 1; return gc_const(av, v);
     514             : }
     515             : 
     516             : /* H[i] = i * (H_i - ((p-1)! + 1) / p) mod p */
     517             : static GEN
     518       33931 : get_GH(GEN F, ulong p)
     519             : {
     520       33931 :   long i, g, ph = p >> 1;
     521       33931 :   GEN H = cgetg(p, t_VECSMALL), v = Flv_inv_p2(p);
     522       33934 :   H[p-1] = (F[p-1] + 1) / p; g = p - H[p-1];
     523    16420541 :   for (i = 1; i <= ph; i++)
     524             :   {
     525    16386619 :     long c = g = Fl_add(g, v[i], p);
     526    16394752 :     H[i] = c = Fl_mul(c, i, p);
     527    16377052 :     H[p-1 - i] = Fl_neg(Fl_add(c, g, p), p);
     528             :   }
     529       33922 :   return H;
     530             : }
     531             : /* p odd, GPB[m+1] = Gamma_p(m/(p-1)) mod p^2 */
     532             : static GEN
     533       33933 : doprecompmodp2(ulong p, ulong p2)
     534             : {
     535       33933 :   GEN F = get_GF(p, p2), H = get_GH(F, p), GPV = cgetg(p, t_VECSMALL);
     536             :   ulong r;
     537    33022669 :   for (r = 1; r < p; r++)
     538             :   {
     539    32988745 :     long t = Fl_mul(F[r], 1 + p * H[r], p2);
     540    32984566 :     GPV[p - r] = odd(r)? t: p2 - t;
     541             :   }
     542       33924 :   return GPV;
     543             : }
     544             : 
     545             : /***************************************************************/
     546             : /*        PART III: Motive initializations depending on p      */
     547             : /***************************************************************/
     548             : static GEN
     549       37597 : get_L1(GEN hgm, long PFM1, long f)
     550             : {
     551       37597 :   GEN VBEDEN = hgm_get_VBEDEN(hgm);
     552       37597 :   GEN VBENUM = hgm_get_VBENUM(hgm), v;
     553       37597 :   long DEG = hgm_get_DEG(hgm), j;
     554       37597 :   v = const_vecsmall(PFM1, f * hgm_get_TT(hgm));
     555      164024 :   for (j = 1; j <= DEG; j++)
     556      126427 :     if (PFM1 % VBEDEN[j] == 0)
     557             :     {
     558      118307 :       long t = ((-(PFM1 / VBEDEN[j]) * VBENUM[j]) % PFM1) + 1;
     559      118307 :       if (t <= 0) t += PFM1;
     560      118307 :       v[t] -= f;
     561             :     }
     562       37597 :   return v;
     563             : }
     564             : /* Stickelberger valuation: L0(p, f, m), power of (-p)
     565             :  * = f * OFFMPOL - sum_{e < f, N < l} gamma_N [N p^e m / (p^f-1)]*/
     566             : /* 1) amortized (all m < p^f-1), good when p >> N log N. Assume f=1 */
     567             : static GEN
     568       33968 : get_L0(GEN hgm, ulong PM1)
     569             : {
     570       33968 :   GEN perm, vt, vc, VL0, VPOLGA = hgm_get_VPOLGA(hgm);
     571       33968 :   long w, S, c, d, j, N, m, l = lg(VPOLGA), D = (l * (l-1)) >> 1;
     572       33968 :   vt = cgetg(D+1, t_VECSMALL);
     573       33967 :   vc = cgetg(D+1, t_VECSMALL);
     574      109098 :   for (N = 2, c = 1; N < l; N++)
     575       75132 :     if (VPOLGA[N])
     576             :     {
     577             :       ulong Q;
     578             :       long q;
     579      195515 :       for (q = 1, Q = 0; q <= N; q++, Q += PM1, c++)
     580             :       {
     581      147137 :         vt[c] = ceildivuu(Q, N);
     582      147135 :         vc[c] = VPOLGA[N];
     583             :       }
     584             :     }
     585       33966 :   setlg(vt,c); setlg(vc,c); perm = vecsmall_indexsort(vt);
     586       33966 :   vt = vecsmallpermute(vt, perm);
     587       33966 :   vc = vecsmallpermute(vc, perm);
     588       33969 :   w = vt[1];
     589      147146 :   for (j = 2, d = 1; j < c; j++)
     590      113177 :     if (vt[j] == w) vc[d] += vc[j];
     591             :     else
     592             :     {
     593       79293 :       d++;
     594       79293 :       vt[d] = w = vt[j];
     595       79293 :       vc[d] = vc[j];
     596             :     }
     597       33969 :   d++; vt[d] = PM1; vc[d] = 0; /* sentinels */
     598       33969 :   VL0 = cgetg(PM1+1, t_VECSMALL); S = hgm_get_OFFMPOL(hgm);
     599      147240 :   for (j = 1; j < d; j++)
     600             :   {
     601    34880338 :     for (m = vt[j]; m < vt[j+1]; m++) VL0[m+1] = S;
     602      113269 :     S -= vc[j+1];
     603             :   }
     604       33971 :   return VL0;
     605             : }
     606             : /* 2) direct computation (single m), good when p << N log N */
     607             : static long
     608      614872 : L0(GEN hgm, ulong p, long f, long PFM1, long m)
     609             : {
     610      614872 :   GEN VPOLGA = hgm_get_VPOLGA(hgm);
     611      614871 :   long e, S = hgm_get_OFFMPOL(hgm) * f, l = lg(VPOLGA);
     612             :   ulong pem;
     613     1858446 :   for (e = 0, pem = m; e < f; e++, pem *= p)
     614             :   {
     615     1243576 :     long N, s = 0;
     616             :     ulong Npem;
     617     9691555 :     for (N = 1, Npem = pem; N < l; N++, Npem += pem)
     618     8447979 :       if (VPOLGA[N]) s += VPOLGA[N] * (long)(Npem / PFM1);
     619     1243576 :     S -= s;
     620             :   }
     621      614870 :   return S;
     622             : }
     623             : 
     624             : static GEN
     625        2947 : get_teich1(GEN VPOLGA, GEN ZP, ulong p)
     626             : {
     627        2947 :   long l = lg(VPOLGA), N;
     628        2947 :   GEN v = zerovec(l - 1);
     629       12712 :   for (N = 2; N < l; N++)
     630        9765 :     if (VPOLGA[N])
     631             :     {
     632        3276 :       long N0; (void)z_lvalrem(N, p, &N0);
     633        3276 :       gel(v, N) = teich(gaddsg(N0, ZP));
     634             :     }
     635        2947 :   return v;
     636             : }
     637             : static GEN
     638        3661 : get_teich(GEN VPOLGA, GEN ZP, ulong p, long f, long PFM1)
     639             : {
     640             :   GEN v, Q;
     641             :   long l, N;
     642        3661 :   if (f == 1) return get_teich1(VPOLGA, ZP, p);
     643         714 :   l = lg(VPOLGA); v = zerovec(l - 1); Q = utoipos(PFM1 / (p - 1));
     644        4277 :   for (N = 2; N < l; N++)
     645        3563 :     if (VPOLGA[N])
     646             :     {
     647        1421 :       long N0; (void)z_lvalrem(N, p, &N0);
     648        1421 :       gel(v,N) = Qp_sqrtn(gdivsg(N0, teich(gaddsg(N0,ZP))), Q, NULL);
     649             :     }
     650         714 :   return v;
     651             : }
     652             : 
     653             : /* compute N^{-[N*s-1-(N*s-1)\p]}, N*s=a/(p^f-1) with a=N*m*p^e */
     654             : static GEN
     655      984431 : gapnpow(GEN T, long p, long f, long PFM1, long N, ulong Nmpe)
     656             : {
     657             :   GEN Vr;
     658             :   long Nm, i, S;
     659             : 
     660      984431 :   if (f == 1) return gpowgs(T, Nmpe);
     661      946295 :   Nm = Fl_neg(Nmpe, PFM1);
     662      946295 :   Vr = cgetg(f + 1, t_VECSMALL);
     663     3372390 :   for (i = 1; i <= f; i++) { Vr[i] = Nm % p; Nm /= p; }
     664     2426095 :   S = Vr[1]; for (i = 1; i < f; i++) S = Vr[f + 1 - i] + p * S;
     665      946295 :   return gdiv(gpowgs(T, S), powuu(N, Vr[1]));
     666             : }
     667             : 
     668             : /* prod_{j = 1}^DEG prod_{e = 0}^{f - 1}
     669             :  *   gamma_p({p^e * (VAL[j] + m/(p^f-1))}) /
     670             :  *   gamma_p({p^e * (VBE[j] + m/(p^f-1))}), a p-unit. 0 <= m < p^f-1 */
     671             : static GEN
     672      235081 : hgmC(GEN VPOLGA, GEN GPV, GEN TEICH, ulong p, long f, ulong PFM1, ulong m, long dfp)
     673             : {
     674      235081 :   GEN r = gen_1;
     675      235081 :   long c, e, N, l = lg(VPOLGA);
     676             :   ulong Nmpe, mpe;
     677      700184 :   for (e = 0, mpe = m; e < f; e++, mpe = Fl_mul(mpe, p, PFM1))
     678     3819344 :     for (N = 1, Nmpe = mpe; N < l; N++, Nmpe = Fl_add(Nmpe, mpe, PFM1))
     679     3354241 :       if ((c = VPOLGA[N]))
     680             :       { /* Gamma_p(frac(N*m*p^e/(p^f-1)))  */
     681     1449522 :         GEN z = utoi(GPV[Nmpe + 1]);
     682     1449508 :         if (N > 1) z = gmul(z, gapnpow(gel(TEICH, N), p, f, PFM1, N, Nmpe));
     683             :         /* z = prod_{0 <= j < N} Gamma_p( frac(p^e * (m/(p^f-1) + j/N)) ) */
     684     1449522 :         r = gmul(r, gpowgs(z, c));
     685             :       }
     686      235073 :   if (typ(r) != t_PADIC) r = cvtop(r, utoipos(p), dfp);
     687      235073 :   return r;
     688             : }
     689             : 
     690             : /* Same modulo p^2, Wm[N] = teich(N^(sign(VPOLGA[N]) * m)) */
     691             : static ulong
     692    16539519 : hgmCmodp2(GEN VPOLGA, GEN Wm, GEN GPV, ulong PFM1, ulong m, ulong p2)
     693             : {
     694    16539519 :   long l = lg(VPOLGA), c, N;
     695    16539519 :   ulong res = 1, Nm;
     696    65307443 :   for (N = 1, Nm = m; N < l; N++, Nm = Fl_add(Nm, m, PFM1))
     697    48836278 :     if ((c = VPOLGA[N]))
     698             :     {
     699             :       ulong z;
     700    37341564 :       if (c > 0) z = GPV[Nm + 1];
     701             :       else
     702             :       {
     703    19171204 :         c = -c;
     704    19171204 :         if (!Nm) z = 1;
     705             :         else
     706             :         {
     707    19151798 :           z = GPV[PFM1 + 1 - Nm];
     708    19151798 :           if (!odd(Nm)) z = p2 - z; /* GPV[Nm+1]^(-1) by reflection formula */
     709             :         }
     710             :       }
     711    38022636 :       if (N > 1) z = Fl_mul(z, Wm[N], p2);
     712    37954105 :       z = Fl_powu(z, c, p2);
     713             :       /* z = (GPV[Nm+1] * teich(N^m))^VPOLGA[N] */
     714    37465781 :       res = res == 1? z: Fl_mul(res, z, p2);
     715             :     }
     716    16413169 :   return res;
     717             : }
     718             : 
     719             : /***************************************************************/
     720             : /*        PART IV: Motive functions depending on p, t          */
     721             : /***************************************************************/
     722             : static long
     723       37896 : get_dfp(GEN hgm, long p, long f)
     724             : {
     725       37896 :   long DEG = hgm_get_DEG(hgm), WT = hgm_get_WT(hgm);
     726       37895 :   long dfp = ceil(log((double)2*DEG) / log((double)p) + f * (WT * 0.5 ) + 1e-5);
     727       37895 :   if (p == 2) dfp++; /* FIXME: why ? */
     728       37895 :   return dfp;
     729             : }
     730             : /* V being a vecsmall, compute all p^TT*hgmC(m)/hgmC(0) for indices in V */
     731             : static GEN
     732        3668 : hgmCall(GEN hgm, ulong p, long f, long dfp, GEN V)
     733             : {
     734        3668 :   GEN v, c, GPV, VL0, VL1, TEICH, ZP = zeropadic_shallow(utoipos(p), dfp);
     735        3668 :   GEN VPOLGA = hgm_get_VPOLGA(hgm);
     736        3668 :   ulong i, PFM1 = upowuu(p, f) - 1, lV = V? lg(V): PFM1+1;
     737        3668 :   long l0, fTT = f * hgm_get_TT(hgm);
     738             : 
     739        3668 :   GPV = doprecomp(p, f, dfp);
     740        3661 :   VL0 = (V && f == 1)? get_L0(hgm, PFM1): NULL;
     741        3661 :   VL1 = get_L1(hgm, PFM1, f);
     742        3661 :   TEICH = get_teich(VPOLGA, ZP, p, f, PFM1);
     743        3661 :   l0 = hgm_get_OFFMPOL(hgm) * f;
     744        3661 :   if (V)
     745             :   {
     746          42 :     v = cgetg(lV, t_VEC);
     747          42 :     i = 1;
     748             :   }
     749             :   else
     750             :   {
     751        3619 :     v = cgetg(lV+1, t_POL); v[1] = evalsigne(1)|evalvarn(0);
     752        3619 :     gel(v,2) = c = powuu(p, fTT); /* m = 0 */
     753        3619 :     i = 2;
     754             :   }
     755      618608 :   for (; i < lV; i++)
     756             :   {
     757      614948 :     long m = V? V[i]: i-1;
     758      614948 :     long L = VL0? VL0[m+1]: L0(hgm, p, f, PFM1, m);
     759      614947 :     long e = L + VL1[m+1];
     760      614947 :     if (!V && e >= dfp) c = gen_0;
     761             :     else
     762             :     { /* FIXME: dfp is fishy in Jordantame, don't trust if V = NULL */
     763      235078 :       pari_sp av = avma;
     764      235078 :       c = hgmC(VPOLGA, GPV, TEICH, p, f, PFM1, m, dfp);
     765      235073 :       setvalp(c, e); if (odd(L ^ l0)) c = gneg(c);
     766      235073 :       c = gerepileupto(av, padic_to_Q(c));
     767             :     }
     768      614947 :     gel(v, V? i: i+1) = c;
     769             :   }
     770        3660 :   return v;
     771             : }
     772             : /* Same mod p^2, f = 1 */
     773             : static GEN
     774       33940 : hgmCallmodp2(GEN hgm, ulong p)
     775             : {
     776       33940 :   GEN C, GPV, VL0, VL1, W1, Wm, VPOLGA = hgm_get_VPOLGA(hgm);
     777       33940 :   long l = lg(VPOLGA), TT = hgm_get_TT(hgm);
     778       33940 :   ulong m, PFM1 = p - 1, p2 = p * p;
     779             : 
     780       33940 :   if (p & HIGHMASK) pari_err_OVERFLOW("hgmCallmodp2");
     781       33933 :   VL0 = get_L0(hgm, PFM1);
     782       33936 :   VL1 = get_L1(hgm, PFM1, 1);
     783       33936 :   W1 = teichmodp2(VPOLGA, p, p2);
     784       33933 :   Wm = shallowcopy(W1);
     785       33933 :   GPV = doprecompmodp2(p, p2);
     786       33924 :   C = cgetg(PFM1+2, t_VECSMALL); C[1] = evalvarn(0);
     787       33930 :   C[2] = TT > 1? 0: (TT == 1? p : 1); /* m = 0 */
     788    32335447 :   for (m = 1; m < PFM1; m++)
     789             :   {
     790    32850799 :     long e = VL0[m + 1] + VL1[m + 1], j;
     791             :     ulong c;
     792    32850799 :     if (e >= 2) c = 0;
     793             :     else
     794             :     {
     795    16447719 :       c = hgmCmodp2(VPOLGA, Wm, GPV, PFM1, m, p2);
     796    16421433 :       if (odd(VL0[m + 1] ^ VL0[1])) c = Fl_neg(c, p2);
     797    15946300 :       if (e == 1) c = (c % p) * p;
     798             :     }
     799    32349380 :     C[m + 2] = c;
     800    96201002 :     for (j = 2; j < l; j++)
     801    63899485 :       if (VPOLGA[j]) Wm[j] = Fl_mul(Wm[j], W1[j], p2);
     802             :   }
     803       29088 :   return C;
     804             : }
     805             : 
     806             : /* 1 / (1-q) ~ 1 + q + ... + q^n; n >= 1 */
     807             : static ulong
     808        3906 : inv(ulong q, long n)
     809             : {
     810        3906 :   ulong z = q + 1;
     811             :   long i;
     812        8064 :   for (i = 1; i < n; i++) z = z * q + 1;
     813        3906 :   return z;
     814             : }
     815             : 
     816             : /* General H function: C(teich(f + O(p^dfp))) / (1 - p^f) */
     817             : static GEN
     818        3906 : hgmH(GEN C, long p, long f, long dfp, GEN t)
     819             : {
     820        3906 :   GEN q = powuu(p, dfp), z;
     821             :   long n;
     822        3906 :   (void)Q_lvalrem(t, p, &t);
     823        3906 :   z = Rg_to_Fp(t, q);
     824        3906 :   z = Zp_teichmuller(z, utoipos(p), dfp, q);
     825        3906 :   z = FpX_eval(C, z, q);
     826        3906 :   n = dfp / f; if (!(dfp % f)) n--;
     827        3906 :   z = Fp_mulu(z, inv(upowuu(p,f), n), q);
     828        3906 :   return Fp_center(z, q,  shifti(q,-1));
     829             : }
     830             : static GEN
     831       33930 : hgmHmodp2(GEN C, ulong p, GEN t)
     832             : {
     833       33930 :   ulong p2 = p * p, wt = Fl_powu(Rg_to_Fl(t, p2), p, p2);
     834       33936 :   return stoi( Fl_center(Fl_mul(Flx_eval(C, wt, p2), 1+p, p2), p2, p2 >> 1) );
     835             : }
     836             : 
     837             : enum { C_OK = 0, C_FAKE, C_BAD, C_TAME0, C_TAME1};
     838             : static GEN hgmU(GEN hgm, long p, long f, GEN t, long dfp);
     839             : static long hgmclass(GEN hgm, long p, GEN t);
     840             : static GEN
     841       37854 : hgmtrace(GEN hgm, long p, long f, GEN t, long c)
     842             : {
     843       37854 :   long dfp = get_dfp(hgm, p, f);
     844       37853 :   if (c == C_FAKE) return hgmU(hgm, p, f, t, dfp);
     845       37566 :   if (f == 1 && dfp <= 2) return hgmHmodp2(hgmCallmodp2(hgm, p), p, t);
     846        3626 :   return hgmH(hgmCall(hgm, p, f, dfp, NULL), p, f, dfp, t);
     847             : }
     848             : 
     849             : static GEN
     850         700 : F2v_factorback(GEN E)
     851             : {
     852         700 :   long i, l = lg(E);
     853         700 :   GEN c = gen_1;
     854        4235 :   for (i = 1; i < l; i++)
     855        3535 :     if (odd(E[i])) c = muliu(c, i);
     856         700 :   return c;
     857             : }
     858             : 
     859             : static long
     860       28420 : Q_krois(GEN T, long p) { return krouu(Rg_to_Fl(T, (p == 2)? 8: p), p); }
     861             : static long
     862       33936 : hgmsign(GEN hgm, long p, GEN t)
     863             : {
     864             :   GEN u;
     865       33936 :   if (odd(hgm_get_WT(hgm))) return 1;
     866       28420 :   t = ginv(t); u = gmul(gsubsg(1, t), hgm_get_SIGNPAR(hgm));
     867       28420 :   return odd(hgm_get_DEG(hgm))? -Q_krois(u, p): Q_krois(gmul(gneg(t), u), p);
     868             : }
     869             : /* conductor of the central character */
     870             : static GEN
     871         259 : hgmcharcond(GEN hgm, GEN t)
     872             : {
     873             :   GEN u;
     874         259 :   if (odd(hgm_get_WT(hgm))) return gen_1;
     875         210 :   t = ginv(t); u = gmul(gsubsg(1, t), hgm_get_SIGNPAR(hgm));
     876         210 :   if (!odd(hgm_get_DEG(hgm))) u = gmul(gneg(t), u);
     877         210 :   if (typ(u) == t_FRAC) u = mulii(gel(u,1), gel(u,2));
     878         210 :   return gequal0(u) ? gen_1 : coredisc(u);
     879             : }
     880             : /* valuations of central character conductor at BAD */
     881             : static GEN
     882          98 : get_achi(GEN hgm, GEN t, GEN BAD)
     883             : {
     884          98 :   long i, l = lg(BAD);
     885          98 :   GEN Nchi = hgmcharcond(hgm, t), v = cgetg(l, t_VECSMALL);
     886         231 :   for (i = 1; i < l; i++) v[i] = Z_lval(Nchi, BAD[i]);
     887          98 :   return v;
     888             : }
     889             : 
     890             : /* [a,a*r, ..., a*r^n] */
     891             : static GEN
     892         420 : upowers_u(ulong r, long n, ulong a)
     893             : {
     894         420 :   long i, l = n + 2, t = a;
     895         420 :   GEN v = cgetg(l, t_VECSMALL);
     896        2562 :   for (i = 1; i < l; i++) { v[i] = t; t *= r; }
     897         420 :   return v;
     898             : }
     899             : static GEN
     900       36876 : powers_u(ulong r, long n, ulong a)
     901             : {
     902       36876 :   long i, l = n + 2;
     903       36876 :   GEN v = cgetg(l, t_VEC), t = utoi(a);
     904      106594 :   for (i = 1; i < l; i++) { gel(v,i) = t; t = muliu(t, r); }
     905       36875 :   return v;
     906             : }
     907             : static GEN
     908       36876 : mkpowers(long p, long d, long WT)
     909             : {
     910             :   ulong q, r;
     911       36876 :   if (WT == 0) q = r = 1;
     912       35700 :   else if (!odd(d)) q = r = upowuu(p, WT);
     913       21889 :   else if (WT == 1) { q = 1; r = p; }
     914       21889 :   else { q = upowuu(p, WT >> 1); r = q*q; if (odd(WT)) r *= p; }
     915       36876 :   return powers_u(r, (d-1)>>1, q);
     916             : }
     917             : 
     918             : /* complete local factor E (of degree d) at p using functional equation
     919             :  * E(T) = SIGN p^(WT*d)/2 T^d E(1/(p^WT*T)) */
     920             : static GEN
     921       36876 : Efuneq(GEN E, long p, long d, long WT, long SIGN, long B)
     922             : {
     923       36876 :   long j = maxss(d - B, 0), k = d + 1 - j, nE = lg(E)-1, l = (d + 1) >> 1;
     924       36876 :   GEN T = cgetg(k + 1, t_VEC), vp = mkpowers(p, d, WT);
     925       39710 :   for (; j < l; j++, k--)
     926             :   {
     927        2835 :     GEN q = gel(vp,l-j);
     928        2835 :     if (SIGN < 0) q = negi(q); /* q = SIGN * p^(WT(d-2*j) / 2) */
     929        2835 :     gel(T, k) = gmul(q, RgX_coeff(E, j));
     930             :   }
     931       39500 :   for (; k >= nE; k--) gel(T, k) = gen_0;
     932      108609 :   for (; k > 0; k--) gel(T, k) = gel(E, k+1);
     933       36875 :   return RgV_to_RgX(T,0);
     934             : }
     935             : 
     936             : static long
     937       37643 : hgmclass(GEN hgm, long p, GEN t)
     938             : {
     939       37643 :   GEN BAD = hgm_get_BAD(hgm);
     940             :   long ap, bp;
     941       37643 :   if (!umodiu(BAD, p))
     942             :   {
     943         679 :     long v = Q_lval(t, p);
     944         679 :     if (v && v + Q_lval(hgm_get_MVALUE(hgm), p) == 0)
     945             :     {
     946         161 :       GEN N = hgmcharcond(hgm, t);
     947         161 :       if (umodiu(N, p)) return C_FAKE; /* a(chi) = 0 */
     948             :     }
     949         553 :     return C_BAD;
     950             :   }
     951       36964 :   if (typ(t) == t_INT)
     952             :   {
     953       23994 :     ap = umodiu(t, p); if (!ap) return C_TAME0;
     954       23966 :     bp = 1;
     955             :   }
     956             :   else
     957             :   {
     958       12970 :     ap = umodiu(gel(t,1), p); if (!ap) return C_TAME0;
     959       12963 :     bp = umodiu(gel(t,2), p); if (!bp) return C_TAME0;
     960             :   }
     961       36866 :   return ap == bp? C_TAME1 : C_OK;
     962             : }
     963             : 
     964             : /* p good or Tame1; return local factor at p: 1/E + O(x^(B+1)); t a t_VEC,
     965             :  * C a t_VECSMALL giving their class. */
     966             : static GEN
     967       36964 : frobpoltrunc(GEN hgm, GEN t, long c, long p, long B, long *pF)
     968             : {
     969       36964 :   long DEG = hgm_get_DEG(hgm), WT = hgm_get_WT(hgm);
     970       36964 :   long D = isint1(t)? (odd(WT) ? DEG - 2 : DEG - 1): DEG;
     971       36964 :   long f, mi, m, q = upowuu(p, WT >> 1);
     972             :   GEN E, s;
     973             : 
     974       36964 :   mi = minss(B, (c == C_FAKE)? D: D >> 1);
     975       36964 :   m = (mi == D && c == C_TAME1 && !odd(DEG))? mi: mi+1;
     976       36964 :   s = cgetg(m+1, t_POL); s[1] = evalsigne(1)|evalvarn(0);
     977       74806 :   for (f = 1; f < m; f++) gel(s, f+1) = negi( hgmtrace(hgm, p, f, t, c) );
     978       36953 :   s = RgX_renormalize_lg(s, m+1);
     979       36952 :   *pF = 0; s = RgXn_expint(s, m);
     980       36953 :   if (mi == D) return s;
     981       36876 :   if (c == C_TAME1)
     982             :   {
     983        2940 :     long SIGN = kroiu(hgm_get_U(hgm), p);
     984        2940 :     if (odd(DEG))
     985         203 :       E = Efuneq(s, p, DEG-1, WT, SIGN, B);
     986             :     else
     987             :     {
     988        2737 :       GEN T = deg1pol_shallow(stoi(- SIGN * q), gen_1, 0);
     989        2737 :       E = RgXn_mul(s, RgXn_inv(T, m), m);
     990        2737 :       E = Efuneq(E, p, DEG - 2, WT, 1, B);
     991        2736 :       if (!gequal1(t) || !odd(WT)) E = gmul(E, T);
     992             :     }
     993        2939 :     *pF = 1;
     994        2939 :     if (!odd(WT))
     995             :     {
     996             :       GEN T, u, t0;
     997         217 :       long v = Q_lvalrem(gsubgs(t, 1), p, &t0);
     998         217 :       if (!odd(v))
     999             :       {
    1000          63 :         if (typ(t0) == t_FRAC) t0 = mulii(gel(t0,1), gel(t0,2));
    1001          63 :         u = coredisc(mulii(t0, hgm_get_U0(hgm)));
    1002          63 :         T = deg1pol_shallow(stoi(-kroiu(u, p) * q), gen_1, 0);
    1003          63 :         E = gmul(E, T); *pF = 0;
    1004             :       }
    1005             :     }
    1006             :   }
    1007             :   else
    1008             :   {
    1009       33936 :     long SIGN = hgmsign(hgm, p, t);
    1010       33936 :     E = Efuneq(s, p, DEG, WT, SIGN, B);
    1011             :   }
    1012       36875 :   return E;
    1013             : }
    1014             : 
    1015             : GEN
    1016          84 : hgmcoef(GEN hgm, GEN t, GEN n)
    1017             : {
    1018          84 :   pari_sp av = avma;
    1019          84 :   GEN T, P, E, F = check_arith_all(n, "hgmcoef");
    1020             :   long i, lP;
    1021             : 
    1022          84 :   if (!checkhgm(hgm)) pari_err_TYPE("hgmcoef", hgm);
    1023          84 :   if (!is_rational_t(typ(t))) pari_err_TYPE("hgmcoef",t);
    1024          77 :   if (hgm_get_SWAP(hgm)) t = ginv(t);
    1025          77 :   if (!F) { F = Z_factor(n); P = gel(F,1); }
    1026             :   else
    1027             :   {
    1028          21 :     P = gel(F,1);
    1029          21 :     if (lg(P) == 1 || signe(gel(P,1)) <= 0) return gen_0;
    1030          21 :     n = typ(n) == t_VEC? gel(n,1): factorback(F);
    1031             :   }
    1032          77 :   if (signe(n) <= 0) pari_err_DOMAIN("hgmcoef", "n", "<=", gen_0, n);
    1033          77 :   E = gel(F,2); lP = lg(P); T = gen_1;
    1034          77 :   T = gen_1;
    1035         161 :   for (i = 1; i < lP; i++)
    1036             :   {
    1037          84 :     long e, p = itos(gel(P, i)), f = itos(gel(E, i)), c = hgmclass(hgm, p, t);
    1038             :     GEN A;
    1039          84 :     if (c == C_BAD) pari_err_IMPL("hgmcoef for bad primes");
    1040          84 :     A = frobpoltrunc(hgm, t, c, p, f, &e);
    1041          84 :     T = gmul(T, RgX_coeff(RgXn_inv(A, f+1), f));
    1042             :   }
    1043          77 :   return gerepilecopy(av, T);
    1044             : }
    1045             : 
    1046             : static GEN
    1047          14 : count2list(GEN E)
    1048             : {
    1049          14 :   long i, j, r, l = lg(E);
    1050          14 :   GEN v = cgetg(zv_sum(E)+1, t_VECSMALL);
    1051          49 :   for (i = j = 1; i < l; i++)
    1052          56 :     if ((r = E[i])) while(r--) v[j++] = i;
    1053          14 :   setlg(v,j); return v;
    1054             : }
    1055             : static GEN hgminit_i(GEN val, GEN vbe);
    1056             : #if 0
    1057             : static long
    1058             : minval(GEN F, long p)
    1059             : {
    1060             :   long i, d = degpol(F), m = LONG_MAX;
    1061             :   for (i = 1; i <= d; i++) m = minss(m, Z_lval(gel(F,i+2), p) / i);
    1062             :   return m;
    1063             : }
    1064             : static GEN
    1065             : cycloE_filter(GEN A, long p)
    1066             : {
    1067             :   long i, j, l = lg(A);
    1068             :   GEN B = shallowcopy(A);
    1069             :   if (p < l)
    1070             :     for (i = p; i < l; i += p)
    1071             :       if (A[i])
    1072             :       {
    1073             :         (void)z_lvalrem(i / p, p, &j);
    1074             :         B[j] += A[i];
    1075             :         B[i] = 0;
    1076             :       }
    1077             :   return B;
    1078             : }
    1079             : /* doesn't work because A/B is not valid hypergeometric data */
    1080             : static GEN
    1081             : eulfacbadnew(GEN hgm, GEN t, long p, long *pe)
    1082             : {
    1083             :   GEN AB = hgm_get_CYCLOE(hgm), gp = utoipos(p);
    1084             :   GEN A = cycloE_filter(gel(AB,1), p);
    1085             :   GEN B = cycloE_filter(gel(AB,2), p);
    1086             :   GEN hgmp = hgminit_i(count2list(A), count2list(B));
    1087             :   GEN cF, E, F = hgmeulerfactor(hgmp, t, p, &E);
    1088             :   long s, v, d = degpol(F), w = hgm_get_WT(hgm);
    1089             :   *pe = 0;
    1090             :   if (!d) return pol_1(0);
    1091             :   cF = leading_coeff(F);
    1092             :   v = Z_lvalrem(cF, p, &cF);
    1093             :   if (!is_pm1(cF)) return pol_1(0);
    1094             :   s = minval(F, p);
    1095             :   if (s) { F = RgX_unscale(F, powis(gp, -s)); v -= s * d; }
    1096             :   if ((2 * v) % d || (!odd(w) && v % d)) return F;
    1097             :   return RgX_unscale(F, powis(gp, w / 2 - v / d));
    1098             : }
    1099             : #endif
    1100             : static GEN Jordantameexpo(GEN hgm, long v, GEN t0, long p, long *pe);
    1101             : static GEN
    1102       37273 : hgmeulerfactorlimit(GEN hgm, GEN t, long p, long d, long flag, long *pe)
    1103             : {
    1104       37273 :   long c = hgmclass(hgm, p, t);
    1105       37272 :   if (c == C_TAME0)
    1106             :   {
    1107          91 :     long v = Q_lvalrem(t, p, &t);
    1108          91 :     return Jordantameexpo(hgm, v, t, p, pe);
    1109             :   }
    1110       37181 :   if (c != C_BAD) return frobpoltrunc(hgm, t, c, p, d, pe);
    1111         301 :   if (flag) { *pe = -1; return gen_0; } else { *pe = 0; return pol_1(0); }
    1112             : }
    1113             : 
    1114             : GEN
    1115         455 : hgmeulerfactor(GEN hgm, GEN t, long p, GEN* pE)
    1116             : {
    1117         455 :   pari_sp av = avma;
    1118             :   long e, B;
    1119             :   GEN P;
    1120         455 :   if (!checkhgm(hgm)) pari_err_TYPE("hgmeulerfactor", hgm);
    1121         455 :   if (!is_rational_t(typ(t))) pari_err_TYPE("hgmeulerfactor",t);
    1122         448 :   if (hgm_get_SWAP(hgm)) t = ginv(t);
    1123         448 :   B = (long)(hgm_get_DEG(hgm) * log(p)) + 1;
    1124         448 :   P = gerepilecopy(av, hgmeulerfactorlimit(hgm, t, p, B, 1, &e));
    1125         434 :   if (pE) *pE = stoi(e);
    1126         434 :   return P;
    1127             : }
    1128             : 
    1129             : /***********************************************************************/
    1130             : /*                        Tame Euler factors                           */
    1131             : /***********************************************************************/
    1132             : /* FIXME: implement properly like RgXn_sqrt */
    1133             : static GEN
    1134          42 : RgXn_sqrtnu(GEN h, long f, long e)
    1135             : {
    1136          42 :   if (f == 1) return h;
    1137           7 :   if (f == 2) return RgXn_sqrt(h, e);
    1138           0 :   return ser2rfrac_i(gsqrtn(RgX_to_ser(h, e + 2), utoipos(f), NULL, 0));
    1139             : }
    1140             : static GEN
    1141          77 : Jordantame(GEN hgm, GEN t0, long m, long p)
    1142             : {
    1143             :   GEN P, T, C, ZP, V;
    1144             :   long d, phim, f, j, c, q, qm, dfp;
    1145             : 
    1146          77 :   if (m == 1)
    1147             :   {
    1148          35 :     long e = hgm_get_WT(hgm) - get_b1(hgm_get_CYCLOE(hgm));
    1149          35 :     return deg1pol_shallow(negi(powuu(p, (e+1) >> 1)), gen_1, 0);
    1150             :   }
    1151          42 :   phim = eulerphiu(m); f = Fl_order(p % m, phim, m);
    1152          42 :   d = phim + 1; /* DEG >= phim >= f */
    1153          42 :   q = upowuu(p, f); qm = (q - 1) / m;
    1154          42 :   V = cgetg(phim + 1, t_VECSMALL);
    1155         133 :   for (j = c = 1; j < m; j++)
    1156          91 :     if (cgcd(j, m) == 1) V[c++] = j * qm;
    1157          42 :   dfp = get_dfp(hgm, p, f);
    1158          42 :   C = hgmCall(hgm, p, f, dfp, V);
    1159          42 :   ZP = zeropadic_shallow(utoipos(p), dfp);
    1160          42 :   T = teich(gadd(t0, ZP)); P = pol_1(0);
    1161         133 :   for (j = 1; j < lg(V); j++)
    1162             :   {
    1163          91 :     GEN Q, c = gmul(gpowgs(T, V[j]), gel(C,j));
    1164          91 :     Q = RgXn_red_shallow(RgX_shift_shallow(RgX_Rg_mul(P, c), f), d);
    1165          91 :     P = RgX_sub(P, Q); /* P * (1 - c x^f) mod X^d */
    1166             :   }
    1167          42 :   return centerlift(RgXn_sqrtnu(P, f, d));
    1168             : }
    1169             : 
    1170             : static GEN
    1171         105 : eulfactameinit(GEN hgm, long v)
    1172             : {
    1173         105 :   GEN C = hgm_get_CYCLOE(hgm);
    1174         105 :   if (!v) pari_err_BUG("hgmeulerfactor [incorrect t in eulfactame]");
    1175         105 :   if (hgm_get_SWAP(hgm)) v = -v;
    1176         105 :   return gel(C, (v < 0)? 2: 1);
    1177             : }
    1178             : static long
    1179          14 : tameexpo(GEN hgm, long v)
    1180             : {
    1181          14 :   GEN W = eulfactameinit(hgm,v);
    1182          14 :   long e, m, l = lg(W);
    1183          42 :   for (m = 1, e = 0; m < l; m++)
    1184          28 :     if (W[m] && v % m == 0) e += eulerphiu(m);
    1185          14 :   return hgm_get_DEG(hgm) - e;
    1186             : }
    1187             : static GEN
    1188          91 : Jordantameexpo(GEN hgm, long v, GEN t0, long p, long *pe)
    1189             : {
    1190          91 :   GEN P = pol_1(0), W = eulfactameinit(hgm, v);
    1191          91 :   long e, m, l = lg(W);
    1192         308 :   for (m = 1, e = 0; m < l; m++)
    1193         217 :     if (W[m] && v % m == 0)
    1194             :     {
    1195          77 :       P = gmul(P, Jordantame(hgm, t0, m, p));
    1196          77 :       e += eulerphiu(m);
    1197             :     }
    1198          91 :   *pe = hgm_get_DEG(hgm) - e; return P;
    1199             : }
    1200             : 
    1201             : /***************************************************************/
    1202             : /*         PART IV.5: Fake wild primes for t                   */
    1203             : /***************************************************************/
    1204             : /* Compute g_q(r)=pi^sq(r)*gauss(r) */
    1205             : static GEN
    1206       57932 : gkgauss(long p, long f, GEN vp, long r, GEN ZP, GEN *sp)
    1207             : {
    1208       57932 :   GEN S = gen_0, P = gen_m1;
    1209       57932 :   long i, qm1 = vp[f+1] - 1;
    1210      272062 :   for (i = 1; i <= f; i++)
    1211             :   {
    1212      214130 :     GEN t = gfrac(sstoQ(r * vp[i], qm1));
    1213      214130 :     S = gadd(S, t);
    1214      214130 :     P = gmul(P, Qp_gamma(gadd(t, ZP)));
    1215             :   }
    1216       57932 :   *sp = gmulsg(p - 1, S); /* t_INT */
    1217       57932 :   return P;
    1218             : }
    1219             : 
    1220             : static GEN
    1221       12292 : hgmG(GEN hgm, long p, long f, GEN vp, long r, GEN ZP)
    1222             : {
    1223       12292 :   GEN S = gen_0, P = gen_1, VPOLGA = hgm_get_VPOLGA(hgm);
    1224             :   long n, c;
    1225       81088 :   for (n = 1; n < lg(VPOLGA); n++)
    1226       68796 :     if ((c = VPOLGA[n]))
    1227             :     {
    1228       57932 :       GEN sq, g = gkgauss(p, f, vp, r*n, ZP, &sq);
    1229       57932 :       S = addii(S, mulsi(c, sq));
    1230       57932 :       P = gmul(P, gpowgs(g, c));
    1231             :     }
    1232             :   /* p - 1 | S */
    1233       12292 :   return gmul(P, powis(utoineg(p), itos(S) / (p - 1)));
    1234             : }
    1235             : 
    1236             : /* multiplicity of -r / (q-1) in beta */
    1237             : static long
    1238       24584 : hgmmulti(GEN B, long q, long r)
    1239             : {
    1240       24584 :   long d = (q-1) / ugcd(r, q-1);
    1241       24584 :   return d >= lg(B)? 0: B[d];
    1242             : }
    1243             : 
    1244             : /* We must have w(M)^r*hgmQ(q,r)=hgmC(q,r)/hgmC(q,0) */
    1245             : static GEN
    1246       12292 : hgmQ(GEN hgm, long p, long f, GEN vp, long r, GEN ZP)
    1247             : {
    1248       12292 :   pari_sp av = avma;
    1249       12292 :   GEN B = gel(hgm_get_CYCLOE(hgm), 2);
    1250       12292 :   long q = vp[f+1], m0 = hgmmulti(B, q, 0), m1 = hgmmulti(B, q, r);
    1251       12292 :   GEN c = powis(utoipos(q), hgm_get_TT(hgm) + m0 - m1);
    1252       12292 :   if (odd(m0)) c = negi(c);
    1253       12292 :   return gerepileupto(av, padic_to_Q(gmul(c, hgmG(hgm, p, f, vp, r, ZP))));
    1254             : }
    1255             : 
    1256             : static GEN
    1257         287 : hgmU(GEN hgm, long p, long f, GEN t, long dfp)
    1258             : {
    1259         287 :   pari_sp av = avma;
    1260         287 :   GEN ZP = zeropadic_shallow(utoipos(p), dfp), vp = upowers_u(p, f, 1);
    1261         287 :   long q = vp[f+1], i;
    1262         287 :   GEN Q = cgetg(q+1, t_POL);
    1263         287 :   Q[1] = evalsigne(1)|evalvarn(0);
    1264       12579 :   for (i = 2; i <= q; i++) gel(Q, i) = hgmQ(hgm, p, f, vp, i-2, ZP);
    1265         287 :   t = p == 2? gen_1: gmul(hgm_get_MVALUE(hgm), t);
    1266         287 :   return gerepileupto(av, hgmH(Q, p, f, dfp, t));
    1267             : }
    1268             : 
    1269             : /***************************************************************/
    1270             : /*         PART V: Utility programs and main init              */
    1271             : /***************************************************************/
    1272             : static GEN
    1273           7 : cycloE2cyclo(GEN A, GEN B)
    1274           7 : { return mkvec2(count2list(A), count2list(B)); }
    1275             : #if 0
    1276             : GEN
    1277             : hgmalphatocyclo(GEN val, GEN vbe)
    1278             : { GEN C = get_CYCLOE(val,vbe); return cycloE2cyclo(gel(C,1), gel(C,2)); }
    1279             : #endif
    1280             : 
    1281             : static GEN
    1282      843850 : allprims(long n, GEN cache)
    1283             : {
    1284             :   long l, i, c;
    1285             :   GEN z, v;
    1286      843850 :   if (gel(cache,n)) return gel(cache,n);
    1287      418635 :   z = coprimes_zv(n); l = lg(z); v = cgetg(l, t_VEC);
    1288     3036782 :   for (i = c = 1; i < l; i++)
    1289     2618147 :     if (z[i]) gel(v, c++) = mkfracss(i, n);
    1290      418635 :   setlg(v, c); return gel(cache,n) = v;
    1291             : }
    1292             : static GEN
    1293      207956 : zv_to_prims(GEN D, GEN cache)
    1294             : {
    1295      207956 :   long i, l = lg(D);
    1296      207956 :   GEN v = cgetg(l, t_VEC);
    1297     1051806 :   for (i = 1; i < l; i++)
    1298             :   {
    1299      843857 :     if (D[i] <= 0) pari_err_TYPE("hgmcyclotoalpha", D);
    1300      843850 :     gel(v, i) = allprims(D[i], cache);
    1301             :   }
    1302      207949 :   return shallowconcat1(v);
    1303             : }
    1304             : static void
    1305      103978 : hgmcyclotoalpha(GEN *pA, GEN *pB)
    1306             : {
    1307      103978 :   GEN v, D = *pA, E = *pB;
    1308      103978 :   if (typ(D) != t_VECSMALL) D = gtovecsmall(D);
    1309      103978 :   if (typ(E) != t_VECSMALL) E = gtovecsmall(E);
    1310      103978 :   v = const_vec(maxss(vecsmall_max(D), vecsmall_max(E)), NULL);
    1311      103978 :   gel(v,1) = mkvec(gen_0);
    1312      103978 :   *pA = zv_to_prims(D, v);
    1313      103978 :   *pB = zv_to_prims(E, v);
    1314      103971 :   if (lg(*pA) != lg(*pB))
    1315           0 :     pari_err_TYPE("hgminit [incorrect lengths]", mkvec2(D,E));
    1316      103971 : }
    1317             : 
    1318             : static GEN
    1319         231 : hgmalphatogamma(GEN val, GEN vbe) { return get_VPOLGA(get_CYCLOE(val, vbe)); }
    1320             : 
    1321             : static void
    1322          21 : hgmgammatocyclo(GEN VPOLGA, GEN *pD, GEN *pE)
    1323             : {
    1324          21 :   long i, cn, cd, l = lg(VPOLGA);
    1325          21 :   GEN d, n, v = zero_zv(l - 1);
    1326          21 :   if (typ(VPOLGA) != t_VECSMALL) VPOLGA = gtovecsmall(VPOLGA);
    1327         112 :   for (i = 1; i < l; i++)
    1328             :   {
    1329          91 :     long e = VPOLGA[i];
    1330          91 :     if (e)
    1331             :     {
    1332          56 :       GEN D = divisorsu(i);
    1333          56 :       long j, lD = lg(D);
    1334         161 :       for (j = 1; j < lD; j++) v[ D[j] ] += e;
    1335             :     }
    1336             :   }
    1337         112 :   for (i = 1, cn = cd = 0; i < l; i++)
    1338             :   {
    1339          91 :     long e = v[i];
    1340          91 :     if (e > 0) cn += e; else cd -= e;
    1341             :   }
    1342          21 :   if (!cn || !cd) pari_err_TYPE("hgmgammatocyclo", VPOLGA);
    1343          14 :   *pE = n = cgetg(cn+1, t_VECSMALL);
    1344          14 :   *pD = d = cgetg(cd+1, t_VECSMALL);
    1345          98 :   for (i = cn = cd = 1; i < l; i++)
    1346             :   {
    1347          84 :     long e = v[i], j;
    1348         189 :     if (e < 0) for (j = 1; j <= -e; j++) d[cd++] = i;
    1349         126 :     else if (e > 0) for (j = 1; j <= e; j++) n[cn++] = i;
    1350             :   }
    1351          14 : }
    1352             : 
    1353             : static void
    1354          21 : hgmgammatoalpha(GEN VPOLGA, GEN *pA, GEN *pB)
    1355          21 : { hgmgammatocyclo(VPOLGA, pA, pB); hgmcyclotoalpha(pA, pB); }
    1356             : 
    1357             : /* A and B sorted */
    1358             : static long
    1359      554932 : zv_intersect(GEN A, GEN B)
    1360             : {
    1361      554932 :   long a = 1, b = 1, lA = lg(A), lB = lg(B);
    1362     1801464 :   while (a < lA && b < lB)
    1363             :   {
    1364     1697794 :     if      (A[a] < B[b]) a++;
    1365     1379231 :     else if (A[a] > B[b]) b++; else return 1;
    1366             :   }
    1367      103670 :   return 0;
    1368             : }
    1369             : static void
    1370         231 : remove_intersect(GEN *pA, GEN *pB)
    1371             : {
    1372         231 :   GEN V, W, A = *pA, B = *pB;
    1373         231 :   long a = 1, b = 1, v = 1, w = 1, lA, lB;
    1374         231 :   *pA = V = cgetg_copy(A, &lA);
    1375         231 :   *pB = W = cgetg_copy(B, &lB);
    1376        1267 :   while (a < lA && b < lB)
    1377             :   {
    1378        1036 :     int s = gcmp(gel(A,a), gel(B,b));
    1379        1036 :     if   (s < 0) gel(V,v++) = gel(A,a++);
    1380         707 :     else if (s > 0) gel(W,w++) = gel(B,b++);
    1381         273 :     else { a++; b++; }
    1382             :   }
    1383         364 :   while (a < lA) gel(V,v++) = gel(A,a++);
    1384         259 :   while (b < lB) gel(W,w++) = gel(B,b++);
    1385         231 :   setlg(V,v); setlg(W,w);
    1386         231 : }
    1387             : 
    1388             : /* al, be normalized, sorted, unknown intersection */
    1389             : static GEN
    1390         231 : albe2u(GEN al, GEN be)
    1391             : {
    1392             :   GEN ga;
    1393         231 :   remove_intersect(&al, &be);
    1394         231 :   ga = hgmalphatogamma(al, be);
    1395         231 :   return F2v_factorback(ga);
    1396             : }
    1397             : static GEN
    1398         469 : get_u(GEN al, GEN be, GEN CYCLOE, GEN VPOLGA, long DEG, long WT)
    1399             : {
    1400         469 :   GEN u, u0 = F2v_factorback(VPOLGA);
    1401         469 :   long b1 = get_b1(CYCLOE);
    1402         469 :   if (odd(DEG))
    1403             :   {
    1404         147 :     be = QV_normalize(RgV_addhalf(be));
    1405         147 :     u = albe2u(al, be);
    1406         147 :     if ((DEG & 3) == 3) u = negi(u);
    1407             :   }
    1408         322 :   else if (odd(WT)) { u = u0; if ((b1 & 3) == 2) u = negi(u); }
    1409             :   else
    1410             :   {
    1411          84 :     al = QV_normalize(RgV_addhalf(al));
    1412          84 :     u = shifti(albe2u(al, be), 1);
    1413          84 :     if (((DEG + b1) & 3) == 1) u = negi(u);
    1414             :   }
    1415         469 :   u0 = shifti(u0, 1); if ((b1 & 3) == 3) u0 = negi(u0);
    1416         469 :   return mkvec2(coredisc(u), u0);
    1417             : }
    1418             : 
    1419             : static long
    1420          35 : zv_sumeuler(GEN v)
    1421             : {
    1422          35 :   long i, l = lg(v);
    1423          35 :   GEN s = gen_0;
    1424         133 :   for (i = 1; i < l; i++)
    1425             :   {
    1426          98 :     if (v[i] <= 0) pari_err_TYPE("hgminit", v);
    1427          98 :     s = addiu(s, eulerphiu(v[i]));
    1428             :   }
    1429          35 :   return itou(s);
    1430             : }
    1431             : 
    1432             : /* (a, b): format (1) if denominator, format (2) if no denominator,
    1433             :  * format (3) if b not vector. */
    1434             : static GEN
    1435         518 : hgminit_i(GEN a, GEN b)
    1436             : {
    1437         518 :   long ta = typ(a), l = lg(a);
    1438         518 :   if (ta != t_VEC && ta != t_VECSMALL) pari_err_TYPE("hgminit", a);
    1439         511 :   if (!b)
    1440             :   {
    1441          84 :     if (l == 1) initab(a, a); /* error */
    1442          77 :     if (ta == t_VECSMALL || RgV_is_ZV(a))
    1443          49 :     {
    1444             :       long i;
    1445          56 :       if (ta != t_VECSMALL) a = vec_to_vecsmall(a);
    1446         161 :       for (i = 1; i < l; i++)
    1447         126 :         if (a[i] <= 0) break;
    1448          56 :       if (i != l)
    1449          21 :         hgmgammatoalpha(a, &a, &b); /* gamma */
    1450             :       else
    1451             :       { /* cyclo */
    1452          35 :         b = const_vecsmall(zv_sumeuler(a), 1);
    1453          35 :         hgmcyclotoalpha(&a, &b);
    1454             :       }
    1455             :     }
    1456             :     else /* alpha */
    1457          21 :       b = zerovec(l - 1);
    1458             :   }
    1459             :   else
    1460             :   {
    1461         427 :     if (typ(b) != ta) pari_err_TYPE("hgminit", b);
    1462         420 :     if (l > 1 && (ta == t_VECSMALL || (RgV_is_ZV(a) && RgV_is_ZV(b))))
    1463         259 :       hgmcyclotoalpha(&a, &b); /* cyclo */
    1464             :   }
    1465         483 :   return initab(a, b);
    1466             : }
    1467             : GEN
    1468         518 : hgminit(GEN val, GEN vbe)
    1469         518 : { pari_sp av = avma; return gerepilecopy(av, hgminit_i(val, vbe)); }
    1470             : 
    1471             : GEN
    1472          21 : hgmalpha(GEN hgm)
    1473             : {
    1474             :   GEN al, be;
    1475          21 :   if (!checkhgm(hgm)) pari_err_TYPE("hgmalpha", hgm);
    1476          14 :   al = hgm_get_VAL(hgm); be = hgm_get_VBE(hgm);
    1477          14 :   if (hgm_get_SWAP(hgm)) swap(al, be);
    1478          14 :   retmkvec2(gcopy(al), gcopy(be));
    1479             : }
    1480             : GEN
    1481          21 : hgmgamma(GEN hgm)
    1482             : {
    1483          21 :   pari_sp av = avma;
    1484             :   GEN g;
    1485          21 :   if (!checkhgm(hgm)) pari_err_TYPE("hgmgamma", hgm);
    1486          14 :   g = hgm_get_VPOLGA(hgm);
    1487          14 :   if (hgm_get_SWAP(hgm)) g = zv_neg(g);
    1488          14 :   return gerepilecopy(av, g);
    1489             : }
    1490             : GEN
    1491           7 : hgmcyclo(GEN hgm)
    1492             : {
    1493           7 :   pari_sp av = avma;
    1494             :   GEN A, B, C;
    1495           7 :   if (!checkhgm(hgm)) pari_err_TYPE("hgmcyclo", hgm);
    1496           7 :   C = hgm_get_CYCLOE(hgm); A = gel(C,1); B = gel(C,2);
    1497           7 :   if (hgm_get_SWAP(hgm)) swap(A, B);
    1498           7 :   return gerepilecopy(av, cycloE2cyclo(A, B));
    1499             : }
    1500             : 
    1501             : GEN
    1502          21 : hgmtwist(GEN hgm)
    1503             : {
    1504          21 :   pari_sp av = avma;
    1505             :   GEN val, vbe;
    1506          21 :   if (!checkhgm(hgm)) pari_err_TYPE("hgmtwist", hgm);
    1507          21 :   val = hgm_get_VAL(hgm);
    1508          21 :   vbe = hgm_get_VBE(hgm);
    1509          21 :   if (hgm_get_SWAP(hgm)) swap(val, vbe);
    1510          21 :   val = sort(RgV_addhalf(val));
    1511          21 :   vbe = sort(RgV_addhalf(vbe));
    1512          21 :   return gerepilecopy(av, initab(val, vbe));
    1513             : }
    1514             : 
    1515             : GEN
    1516         161 : hgmparams(GEN hgm)
    1517             : {
    1518         161 :   pari_sp av = avma;
    1519             :   GEN H, M;
    1520             :   long TT, DEG, WT;
    1521         161 :   if (!checkhgm(hgm)) pari_err_TYPE("hgmparams", hgm);
    1522         161 :   H = zx_to_ZX(hgm_get_HODGE(hgm));
    1523         161 :   TT = hgm_get_TT(hgm); DEG = hgm_get_DEG(hgm);
    1524         161 :   WT = hgm_get_WT(hgm); M = hgm_get_MVALUE(hgm);
    1525         161 :   return gerepilecopy(av, mkvec4(utoipos(DEG), utoi(WT),
    1526             :                                  mkvec2(H,stoi(TT)), M));
    1527             : }
    1528             : 
    1529             : /* symmetry at 1? */
    1530             : long
    1531           7 : hgmissymmetrical(GEN hgm)
    1532             : {
    1533             :   GEN A, B, C;
    1534             :   long a, i, j, lA, lB;
    1535           7 :   if (!checkhgm(hgm)) pari_err_TYPE("hgmissymmetrical", hgm);
    1536           7 :   C = hgm_get_CYCLOE(hgm);
    1537           7 :   if (odd(hgm_get_DEG(hgm))) return 0;
    1538           7 :   A = gel(C,1); lA = lg(A);
    1539           7 :   B = gel(C,2); lB = lg(B);
    1540          21 :   for (i = 1; i < lA; i++) if ((a = A[i]))
    1541             :   {
    1542           7 :     switch(i & 3)
    1543             :     { /* polcyclo(i)[-X] = polcyclo(j) */
    1544           0 :       case 0: j = i; break;
    1545           7 :       case 2: j = i >> 1; break;
    1546           0 :       default: j = i << 1; break;
    1547             :     }
    1548           7 :     if (j >= lB || B[j] != a) return 0;
    1549             :   }
    1550           7 :   return 1;
    1551             : }
    1552             : 
    1553             : /***************************************************************/
    1554             : /*         PART VI: Experimental euler factors                 */
    1555             : /***************************************************************/
    1556             : /* FIXME: one prime at a time */
    1557             : static GEN
    1558      227731 : hgmmodif(GEN an, GEN PPol)
    1559             : {
    1560      227731 :   pari_sp av = avma;
    1561      227731 :   long i, L = lg(an), lP = lg(PPol);
    1562             : 
    1563      674401 :   for (i = 1; i < lP; i++)
    1564             :   {
    1565      446670 :     GEN E = gel(PPol, i), pol = gel(E, 2);
    1566      446670 :     long d = degpol(pol);
    1567      446670 :     if (d)
    1568             :     {
    1569      433629 :       GEN v = vec_ei(L, 1);
    1570      433629 :       long j, p = itos(gel(E, 1)), pj = p;
    1571     1384271 :       for (j = 1; j <= d && pj <= L; j++, pj *= p)
    1572      950642 :         gel(v, pj) = RgX_coeff(pol, j);
    1573      433629 :       an = dirdiv(an, v);
    1574             :     }
    1575             :   }
    1576      227731 :   return gerepileupto(av, an);
    1577             : }
    1578             : 
    1579             : /***************************************************************/
    1580             : /*             PART VII: Make tables of HGMs                   */
    1581             : /***************************************************************/
    1582             : 
    1583             : static int
    1584         378 : ok_part(GEN v, GEN W)
    1585             : {
    1586         378 :   long l = lg(v), j;
    1587        1617 :   for (j = 1; j < l; j++)
    1588        1435 :     if (!gel(W,v[j])) return 0;
    1589         182 :   return 1;
    1590             : }
    1591             : static GEN
    1592          21 : mkphi()
    1593             : {
    1594          21 :   GEN W = const_vec(20, NULL);
    1595          21 :   gel(W,1) = mkvecsmall2(1,2);
    1596          21 :   gel(W,2) = mkvecsmall3(3,4,6);
    1597          21 :   gel(W,4) = mkvecsmall4(5,8,10,12);
    1598          21 :   gel(W,6) = mkvecsmall4(7,9,14,18);
    1599          21 :   gel(W,8) = mkvecsmall5(15,16,20,24,30);
    1600          21 :   gel(W,10)= mkvecsmall2(11,22);
    1601          21 :   gel(W,12)= mkvecsmalln(6, 13L,21L,26L,28L,36L,42L);
    1602          21 :   gel(W,16)= mkvecsmalln(6, 17L,32L,34L,40L,48L,60L);
    1603          21 :   gel(W,18)= mkvecsmall4(19,27,38,54);
    1604          21 :   gel(W,20)= mkvecsmall5(25,33,44,50,66);
    1605          21 :   return W;
    1606             : }
    1607             : 
    1608             : static GEN
    1609         182 : mkal(GEN p, GEN W)
    1610             : {
    1611             :   GEN res, V;
    1612         182 :   long l = lg(p), lV, i, j;
    1613         182 :   V = cgetg(l, t_VECSMALL);
    1614         987 :   for (i = 1; i < l; i++) V[i] = lg(gel(W, p[i])) - 1;
    1615         182 :   V = cyc2elts(V); lV = lg(V);
    1616         182 :   res = cgetg(1, t_VEC);
    1617       20769 :   for (j = 1; j < lV; j++)
    1618             :   {
    1619       20587 :     GEN t = cgetg(l, t_VECSMALL), v = gel(V,j);
    1620      157185 :     for (i = 1; i < l; i++) t[i] = umael2(W, p[i], v[i]+1);
    1621       20587 :     vecsmall_sort(t); if (!RgV_isin(res, t)) res = vec_append(res, t);
    1622             :   }
    1623         182 :   return res;
    1624             : }
    1625             : static GEN
    1626          21 : mkallal(long n)
    1627             : {
    1628          21 :   GEN W = mkphi(), p = partitions(n, NULL, NULL);
    1629          21 :   long i, c, l = lg(p);
    1630         399 :   for (i = c = 1; i < l; i++)
    1631         378 :     if (ok_part(gel(p,i), W)) gel(p, c++) = mkal(gel(p,i), W);
    1632          21 :   setlg(p, c); return shallowconcat1(p);
    1633             : }
    1634             : 
    1635             : static GEN
    1636          21 : mkalbe(long n)
    1637             : {
    1638          21 :   GEN w, L = mkallal(n);
    1639          21 :   long i, j, c, l = lg(L);
    1640          21 :   w = cgetg(l * (l / 2) + 1, t_VEC);
    1641        3941 :   for (i = c = 1; i < l; i++)
    1642             :   {
    1643        3920 :     GEN A = gel(L, i);
    1644        3920 :     long a = A[lg(A)-1];
    1645      558852 :     for (j = i + 1; j < l; j++)
    1646             :     {
    1647      554932 :       GEN B = gel(L, j);
    1648      554932 :       if (!zv_intersect(A, B))
    1649      179081 :         gel(w,c++) = (A[1]==1 || (B[1]!=1 && a > B[lg(B)-1]))?
    1650      179081 :                      mkvec2(B,A): mkvec2(A,B);
    1651             :     }
    1652             :   }
    1653          21 :   setlg(w,c); return w;
    1654             : }
    1655             : 
    1656             : static long
    1657      103670 : cyclowt(GEN a, GEN b)
    1658             : {
    1659      103670 :   pari_sp av = avma;
    1660             :   long TT;
    1661      103670 :   hgmcyclotoalpha(&a, &b);
    1662      103670 :   return gc_long(av, degpol(hodge(a, b, &TT)));
    1663             : }
    1664             : 
    1665             : GEN
    1666          21 : hgmbydegree(long n)
    1667             : {
    1668          21 :   pari_sp av = avma;
    1669          21 :   GEN w = cgetg(n+1, t_VEC), c = const_vecsmall(n, 1), v = mkalbe(n);
    1670          21 :   long i, l = lg(v);
    1671         154 :   for (i = 1; i <= n; i++) gel(w,i) = cgetg(l,t_VEC);
    1672      103691 :   for (i = 1; i < l; i++)
    1673             :   {
    1674      103670 :     GEN z = gel(v,i);
    1675      103670 :     long k = cyclowt(gel(z,1), gel(z,2)) + 1;
    1676      103670 :     gmael(w, k, c[k]++) = z;
    1677             :   }
    1678         154 :   for (i = 1; i <= n; i++) setlg(gel(w,i), c[i]);
    1679          21 :   return gerepilecopy(av, w);
    1680             : }
    1681             : 
    1682             : /***************************************************************/
    1683             : /*             PART VIII: L-function data                      */
    1684             : /***************************************************************/
    1685             : /* BAD sorted t_VECSMALL */
    1686             : static GEN
    1687         217 : removebad(GEN v, GEN BAD)
    1688             : {
    1689         217 :   long i, c, l = lg(v);
    1690         217 :   GEN w = cgetg(l, t_VECSMALL);
    1691         434 :   for (i = c = 1; i < lg(v); i++)
    1692         217 :     if (!zv_search(BAD, v[i])) w[c++] = v[i];
    1693         217 :   setlg(w, c); return w;
    1694             : }
    1695             : static GEN
    1696         609 : primedivisors(GEN t)
    1697         609 : { GEN fa = absZ_factor(t); return ZV_to_zv(gel(fa,1)); }
    1698             : 
    1699             : static GEN
    1700         385 : gacfac(long p, long m, long c)
    1701             : {
    1702         385 :   long i, l = p + m + 1;
    1703         385 :   GEN v = cgetg(l, t_VECSMALL);
    1704         728 :   for (i = 1; i <= p; i++) v[i] = -c;
    1705         658 :   for (     ; i <  l; i++) v[i] = 1 - c;
    1706         385 :   return v;
    1707             : }
    1708             : 
    1709             : static GEN
    1710         203 : hgmfindvga(GEN hgm, GEN t)
    1711             : {
    1712         203 :   GEN v, HODGE = hgm_get_HODGE(hgm);
    1713         203 :   long WT = degpol(HODGE), WT2 = (WT - 1) >> 1, i, c, h, fl = gequal1(t);
    1714         203 :   v = cgetg(WT + 2, t_VEC);
    1715         413 :   for (i = 0, c = 1; i <= WT2; i++)
    1716             :   {
    1717         210 :     if ((h = HODGE[i+2]))
    1718             :     {
    1719         210 :       if (fl && 2*i == WT - 1) h--;
    1720         210 :       if (h) gel(v, c++) = gacfac(h, h, i);
    1721             :     }
    1722             :   }
    1723         203 :   if (!odd(WT))
    1724             :   {
    1725         182 :     long h = HODGE[WT2+3], hp = h >> 1, hm;
    1726         182 :     if (!fl)
    1727             :     {
    1728         182 :       hm = hp;
    1729         182 :       if (odd(h))
    1730         154 :       { if (gsigne(t) >= 0 && gcmpgs(t, 1) <= 0) hm++; else hp++; }
    1731             :       else
    1732          28 :       { if (gcmpgs(t, 1) > 0) { hp++; hm--; } }
    1733         182 :       if (odd(hgm_get_TT(hgm) + WT2 + 1)) lswap(hp, hm);
    1734             :     }
    1735             :     else
    1736             :     {
    1737           0 :       h--; hm = h - hp;
    1738           0 :       if (odd(h) && odd(hgm_get_TT(hgm) + WT2 + 1)) lswap(hp, hm);
    1739             :     }
    1740         182 :     if (hp || hm) gel(v, c++) = gacfac(hp, hm, WT2 + 1);
    1741             :   }
    1742         203 :   if (c == 1) return cgetg(1, t_VECSMALL);
    1743         203 :   setlg(v, c); v = shallowconcat1(v);
    1744         203 :   vecsmall_sort(v); return v;
    1745             : }
    1746             : 
    1747             : /* Return [VGA, k, BAD, COND] where VGA as in lfun, k non motivic
    1748             :  * weight (s -> k - s), BAD is the list of wild primes and Euler factors,
    1749             :  * COND is the tame part of the conductor */
    1750             : static GEN
    1751         203 : hgmlfuninfty(GEN hgm, GEN t)
    1752             : {
    1753         203 :   pari_sp av = avma;
    1754         203 :   GEN VGA = hgmfindvga(hgm, t), T0, T1, v;
    1755         203 :   GEN COND, t1 = gsubgs(t, 1), BAD = primedivisors(hgm_get_BAD(hgm));
    1756         203 :   long k = hgm_get_WT(hgm) + 1, i, j, l;
    1757             : 
    1758         203 :   if (isintzero(t1)) T1 = cgetg(1, t_VECSMALL);
    1759             :   else
    1760             :   {
    1761         196 :     GEN fa = absZ_factor(numer_i(t1)), P = gel(fa,1), E = gel(fa,2);
    1762         196 :     if (!odd(k)) T1 = removebad(ZV_to_zv(P), BAD);
    1763             :     else
    1764             :     {
    1765         182 :       long j, lP = lg(P);
    1766         182 :       T1 = cgetg(lP, t_VECSMALL);
    1767         350 :       for (i = j = 1; i < lP; i++)
    1768             :       {
    1769         168 :         long p = itos(gel(P,i));
    1770         168 :         if (mpodd(gel(E,i)) && !zv_search(BAD, p)) T1[j++] = p;
    1771             :       }
    1772         182 :       setlg(T1,j);
    1773             :     }
    1774             :   }
    1775         203 :   COND = zv_prod_Z(T1);
    1776         203 :   T0 = removebad(gconcat(primedivisors(numer_i(t)),
    1777             :                          primedivisors(denom_i(t))), BAD);
    1778         217 :   for (i = 1; i < lg(T0); i++)
    1779             :   {
    1780          14 :     long p = T0[i], e = tameexpo(hgm, Q_lval(t, p));
    1781          14 :     COND = mulii(powuu(p, e), COND);
    1782             :   }
    1783         203 :   l = lg(BAD); v = cgetg(l, t_VEC);
    1784         490 :   for (i = j = 1; i < l; i++) /* FIXME: precompute wild Euler factors */
    1785         287 :     if (hgmclass(hgm, BAD[i], t) == C_BAD)
    1786         252 :       gel(v,j++) = mkvec2(utoipos(BAD[i]), pol_1(0));
    1787         203 :   setlg(v, j); return gerepilecopy(av, mkvec4(VGA, utoi(k), v, COND));
    1788             : }
    1789             : 
    1790             : /***************************************************************/
    1791             : /*             PART IX: GUESS OF WILD PRIMES                   */
    1792             : /***************************************************************/
    1793             : static GEN
    1794        1708 : vecmellin(GEN L, GEN K, GEN t, GEN td, GEN vroots, long bitprec)
    1795             : {
    1796        1708 :   long i, N = lfunthetacost(L, t, 0, bitprec);
    1797        1708 :   GEN v = cgetg(N+1, t_VEC);
    1798      427231 :   for (i = 1; i <= N; i++)
    1799      425523 :     gel(v,i) = gammamellininvrt(K, gmul(td, gel(vroots,i)), bitprec);
    1800        1708 :   return v;
    1801             : }
    1802             : 
    1803             : /* List of Weil polynomials for prime p of degree d and weight w */
    1804             : static GEN
    1805        1610 : listweil_i(ulong d, ulong p, ulong w)
    1806             : {
    1807             :   GEN V;
    1808        1610 :   if (d == 0) return mkvec(pol_1(0));
    1809        1225 :   if (odd(d))
    1810             :   {
    1811             :     GEN P;
    1812         644 :     if (odd(w)) return cgetg(1, t_VEC);
    1813         420 :     V = listweil_i(d - 1, p, w);
    1814         420 :     P = monomial(powuu(p, w / 2), 1, 0);
    1815         420 :     return shallowconcat(gmul(gsubsg(1, P), V), gmul(gaddsg(1, P), V));
    1816             :   }
    1817         581 :   if (d == 2)
    1818             :   {
    1819         553 :     long q = upowuu(p, w), s = usqrt(4 * q), i;
    1820         553 :     GEN Q = utoi(q);
    1821         553 :     V = cgetg(2*s + 3, t_VEC);
    1822        4844 :     for (i = 1; i <= 2*s + 1; i++) gel(V,i) = mkpoln(3, Q, stoi(1+s-i), gen_1);
    1823         553 :     gel(V, 2*s + 2) = mkpoln(3, negi(Q), gen_0, gen_1);
    1824         553 :     return V;
    1825             :   }
    1826          28 :   if (d == 4)
    1827             :   {
    1828          28 :     long q = upowuu(p, w), s = usqrt(16 * q), a, c, is2 = usqrt(4 * q);
    1829          28 :     double s2 = 2 * sqrt((double)q);
    1830          28 :     GEN Q2 = sqru(q), W, A, mA, tmp, Q;
    1831          28 :     V = cgetg(s + 3, t_VEC);
    1832         252 :     for (a = 0; a <= s; a++)
    1833             :     {
    1834         224 :       long b, m = ceil(a * s2) - 2 * q, M = ((a * a) >> 2) + 2 * q;
    1835         224 :       GEN AQ = stoi(a * q), mAQ = stoi(-a * q);
    1836         224 :       A = stoi(a); mA = stoi(-a);
    1837         224 :       W = cgetg(2*(M - m) + 3, t_VEC);
    1838        1876 :       for (b = m, c = 1; b <= M; b++)
    1839             :       {
    1840        1652 :         if (a) gel(W, c++) = mkpoln(5, Q2, mAQ, stoi(b), mA, gen_1);
    1841        1652 :         gel(W, c++) = mkpoln(5, Q2, AQ, stoi(b), A, gen_1);
    1842             :       }
    1843         224 :       setlg(W, c); gel(V, a + 1) = W;
    1844             :     }
    1845          28 :     W = cgetg(2 * is2 + 2, t_VEC);
    1846          28 :     tmp = mkpoln(3, stoi(-q), gen_0, gen_1);
    1847          28 :     Q = utoipos(q);
    1848         147 :     for (a = 0, c = 1; a <= is2; a++)
    1849             :     {
    1850         119 :       A = stoi(a); mA = stoi(-a);
    1851         119 :       if (a) gel(W, c++) = gmul(tmp, mkpoln(3, Q, mA, gen_1));
    1852         119 :       gel(W, c++) = gmul(tmp, mkpoln(3, Q, A, gen_1));
    1853             :     }
    1854          28 :     setlg(W, c); gel(V, s + 2) = W;
    1855          28 :     return shallowconcat1(V);
    1856             :   }
    1857           0 :   pari_err_IMPL("d > 5 in listweil");
    1858             :   return NULL; /* LCOV_EXCL_LINE */
    1859             : }
    1860             : 
    1861             : /* Same, weight <= w, by decreasing weight */
    1862             : static GEN
    1863         490 : listweilallw_i(ulong d, ulong p, ulong w)
    1864             : {
    1865         490 :   long i = d? w: 0;
    1866         490 :   GEN V = cgetg(i+2, t_VEC);
    1867        1680 :   for (; i >= 0; i--) gel(V,i+1) = listweil_i(d, p, i);
    1868         490 :   return shallowconcat1(V);
    1869             : }
    1870             : 
    1871             : static long
    1872      127344 : sumdeg(GEN W, GEN M)
    1873             : {
    1874      127344 :   long i, l = lg(M), s = 0;
    1875      375851 :   for (i = 1; i < l; i++) s += degpol(gmael(W,i,M[i]+1));
    1876      127344 :   return s;
    1877             : }
    1878             : 
    1879             : static GEN
    1880          35 : mul(GEN a, GEN b, ulong lim)
    1881             : {
    1882          35 :   long na = lg(a)-1, nb = lg(b)-1, i, j, n;
    1883          35 :   GEN c = cgetg(na * nb + 1, t_VECSMALL);
    1884         378 :   for (i = n = 1; i <= na; i++)
    1885        2597 :     for (j = 1; j <= nb; j++)
    1886             :     {
    1887        2254 :       ulong m = umuluu_or_0(a[i], b[j]);
    1888        2254 :       if (m && m <= lim) c[n++] = m;
    1889             :     }
    1890          35 :   setlg(c, n); return c;
    1891             : }
    1892             : static GEN
    1893          98 : listcond(GEN BAD, GEN achi, ulong min, ulong max)
    1894             : {
    1895          98 :   long i, j, l = lg(BAD);
    1896          98 :   GEN P = cgetg(l, t_VEC), z;
    1897         231 :   for (i = 1; i < l; i++)
    1898             :   {
    1899         133 :     long p = BAD[i], v = achi[i];
    1900         133 :     gel(P,i) = upowers_u(p, ulogint(max, p) - v, upowuu(p, v));
    1901             :   }
    1902         133 :   z = gel(P,1); for (i = 2; i < l; i++) z = mul(z, gel(P,i), max);
    1903          98 :   if (min > 1)
    1904             :   {
    1905           0 :     l = lg(z);
    1906           0 :     for (i = j = 1; i < l; i++)
    1907           0 :       if ((ulong)z[i] >= min) z[j++] = z[i];
    1908           0 :     setlg(z, j);
    1909             :   }
    1910          98 :   vecsmall_sort(z); return z;
    1911             : }
    1912             : 
    1913             : /* Artificial lfundiv by zeta(s - (k - 1) / 2). */
    1914             : static GEN
    1915          35 : lfundivraw(GEN L)
    1916             : {
    1917          35 :   long k = itos(ldata_get_k(L)), i;
    1918          35 :   GEN v = ldata_get_gammavec(L);
    1919          35 :   i = RgV_isin(v, stoi(-(k - 1) >> 2));
    1920          35 :   if (!i) pari_err_BUG("lfundivraw");
    1921          35 :   L = shallowcopy(L);
    1922          35 :   gel(L, 3) = vecsplice(v, i);
    1923          35 :   setlg(L, 7); return L;
    1924             : }
    1925             : 
    1926             : /* Compute forvec vectors from v, sorted by increasing total degree */
    1927             : static GEN
    1928          98 : forvecsort(GEN vF, GEN v)
    1929             : {
    1930          98 :   GEN w, E = cyc2elts(v);
    1931          98 :   long i, l = lg(E);
    1932          98 :   w = cgetg(l, t_VECSMALL);
    1933      127442 :   for (i = 1; i < l; i++) w[i] = sumdeg(vF, gel(E, i));
    1934          98 :   return vecpermute(E, vecsmall_indexsort(w)); /* by increasing total degree */
    1935             : }
    1936             : 
    1937             : static GEN
    1938          98 : lfunhgmwild(GEN L, GEN H, GEN t, GEN BAD, long pole, GEN hint, long bitprec)
    1939             : {
    1940             :   GEN v, K, t0, t0r, t0ir, t0i, t0k, N0, vM, vD, val, PPOL, vF, achi;
    1941          98 :   long d, lM, iN, iM, i, k, k2, prec = nbits2prec(bitprec), lB = lg(BAD);
    1942             :   long BADprod, limdeg;
    1943          98 :   ulong minN = 1, maxN = 2048;
    1944             : 
    1945          98 :   v = cgetg(lB, t_VECSMALL); PPOL = cgetg(lB, t_VEC);
    1946         231 :   for (i = 1; i < lB; i++)
    1947             :   {
    1948         133 :     v[i] = itou( gmael(BAD,i,1) );
    1949         133 :     gel(PPOL,i) = shallowcopy(gel(BAD,i));
    1950             :   }
    1951          98 :   BAD = v;
    1952          98 :   BADprod = zv_prod(BAD);
    1953          98 :   achi = get_achi(H, t, BAD);
    1954          98 :   if (pole) L = lfundivraw(L);
    1955          98 :   limdeg = d = ldata_get_degree(L);
    1956          98 :   N0 = ldata_get_conductor(L);
    1957          98 :   if (hint) switch(typ(hint))
    1958             :   {
    1959             :     long l;
    1960           0 :     case t_INT:
    1961           0 :       limdeg = itos(hint);
    1962           0 :       if (limdeg < 0 || limdeg > d) pari_err_TYPE("lfunhgm [bad hint]", hint);
    1963           0 :       break;
    1964           0 :     case t_VEC:
    1965           0 :       l = lg(hint);
    1966           0 :       if (l > 1 && l < 4)
    1967             :       {
    1968           0 :         GEN t = gel(hint,1), r;
    1969           0 :         if (typ(t) != t_INT || signe(t) <= 0)
    1970           0 :           pari_err_TYPE("lfunhgm [bad hint]", hint);
    1971           0 :         t = dvmdii(t, N0, &r);
    1972           0 :         if (r != gen_0)
    1973           0 :           pari_err_TYPE("lfunhgm [bad hint]", hint);
    1974           0 :         minN = maxN = itou(t);
    1975           0 :         if (l == 3)
    1976             :         {
    1977           0 :           t = gel(hint,2);
    1978           0 :           if (typ(t) != t_INT || signe(t) < 0 || cmpis(t, d) > 0)
    1979           0 :             pari_err_TYPE("lfunhgm [bad hint]", hint);
    1980           0 :           limdeg = itos(t);
    1981             :         }
    1982             :       }
    1983             :   }
    1984          98 :   k = itos(ldata_get_k(L)); k2 = (k-1) >> 1;
    1985          98 :   K = gammamellininvinit(ldata_get_gammavec(L), 0, bitprec + 32);
    1986          98 :   t0 = sstoQ(11, 10); t0i = ginv(t0); t0k = gpowgs(t0, k);
    1987          98 :   t0r = gpow(t0, sstoQ(2,d), prec); t0ir = ginv(t0r);
    1988             :   /* vF[i] list of Weil polynomials F for p = BAD[i],
    1989             :    * F = prod_j Fj, |reciprocal root of Fj|^2 = p^(w + 1 - nj)
    1990             :    * 2v_p(lc(F)) = deg F * (w + 1) - sum_j deg Fj * nj; */
    1991          98 :   vF = cgetg(lB, t_VEC);
    1992          98 :   vD = cgetg(lB, t_VEC); /* vD[k][l] = sum_j deg Fj * nj for F = vF[k][l] */
    1993          98 :   v = cgetg(lB, t_VECSMALL);
    1994         231 :   for (i = 1; i < lB; i++)
    1995             :   {
    1996         133 :     GEN W = cgetg(limdeg+2, t_VEC), D;
    1997         133 :     long p = BAD[i], j, lW;
    1998         623 :     for (j = 0; j <= limdeg; j++) gel(W,j+1) = listweilallw_i(j, p, k-1);
    1999         133 :     gel(vF, i) = W = shallowconcat1(W);
    2000         133 :     lW = lg(W); v[i] = lW-1;
    2001         133 :     gel(vD,i) = D = cgetg(lW, t_VEC);
    2002       10136 :     for (j = 1; j < lW; j++)
    2003             :     {
    2004       10003 :       GEN F = gel(W,j);
    2005       10003 :       D[j] = degpol(F) * k - 2 * Z_lval(leading_coeff(F), p);
    2006             :     }
    2007             :   }
    2008          98 :   vM = forvecsort(vF, v); lM = lg(vM);
    2009          98 :   if (DEBUGLEVEL) { err_printf(" lM = %ld ", lM); err_flush(); }
    2010          98 :   L = shallowcopy(L);
    2011          98 :   val = cgetg(lB, t_VECSMALL);
    2012             :   for(;;)
    2013           0 :   {
    2014             :     GEN z, vroots, an0, vN;
    2015             :     long lN, lim;
    2016             : 
    2017          98 :     z = listcond(BAD, achi, minN, maxN);
    2018          98 :     if (maxN == minN) /* from hint */
    2019             :     {
    2020           0 :       minN = 1; /* in case it fails */
    2021           0 :       maxN--;
    2022             :     }
    2023             :     else
    2024             :     {
    2025          98 :       minN = maxN+1;
    2026          98 :       maxN <<= 2;
    2027             :     }
    2028          98 :     lN = lg(z);
    2029          98 :     if (lN == 1) continue;
    2030          98 :     vN = cgetg(lN, t_VEC);
    2031        1785 :     for (i = 1; i < lN; i++) gel(vN, i) = muliu(N0, z[i]);
    2032          98 :     gel(L, 5) = gel(vN, lN - 1);
    2033          98 :     if (DEBUGLEVEL) err_printf("\nmaxN = %ld\n", itos(gel(L,5)));
    2034          98 :     lim = lfunthetacost(L, t0i, 0, bitprec);
    2035          98 :     vroots = dirpowers(lim, sstoQ(2, d), prec + EXTRAPREC64);
    2036          98 :     an0 = hgmcoefs(H, t, lim);
    2037          98 :     if (pole)
    2038             :     {
    2039          35 :       GEN w = vecpowuu(lim, k2);
    2040       44450 :       for (i = 1; i <= lim; i++)
    2041       44415 :         if (cgcd(i, BADprod) > 1) gel(w, i) = gen_0;
    2042          35 :       an0 = dirdiv(an0, w);
    2043             :     }
    2044         854 :     for (iN = 1; iN < lN; iN++)
    2045             :     {
    2046         854 :       pari_sp av0 = avma, av2;
    2047         854 :       GEN vecK, vecKi, N = gel(vN,iN), isqN = gpow(N, sstoQ(-1,d), prec);
    2048         854 :       if (DEBUGLEVEL) { err_printf(" %ld ", itos(N)); err_flush(); }
    2049         854 :       gel(L,5) = N;
    2050         854 :       vecK = vecmellin(L, K, t0, gmul(t0r,  isqN), vroots, bitprec);
    2051         854 :       vecKi= vecmellin(L, K, t0i,gmul(t0ir, isqN), vroots, bitprec);
    2052        2310 :       for (i = 1; i < lB; i++) val[i] = Z_lval(N, BAD[i]);
    2053         854 :       setlg(an0, lg(vecKi)); av2 = avma;
    2054     2305926 :       for (iM = 1; iM < lM; iM++, set_avma(av2))
    2055             :       {
    2056     2305170 :         GEN M = gel(vM, iM), an, eno, den;
    2057     3214589 :         for (i = 1; i < lB; i++)
    2058             :         {
    2059     2986858 :           GEN F = gmael(vF, i, M[i]+1);
    2060     2986858 :           long D, dF = degpol(F), a = val[i];
    2061     2986858 :           if (a < d - dF) break;
    2062     2797417 :           if (a)
    2063             :           {
    2064     2606261 :             if (d == dF) break;
    2065     1163771 :             D = mael(vD, i, M[i]+1); /* sum_j deg Fj nj */
    2066     1163771 :             if (D == d && d - dF != a) break;
    2067             :             /* a = a(chi) + d - deg F - 1 */
    2068     1100575 :             if (D == d-1 && a != d - dF + achi[i] - 1 && !gequal1(t)) break;
    2069             :           }
    2070      909419 :           gmael(PPOL, i, 2) = F;
    2071             :         }
    2072     2305170 :         if (i < lB) continue;
    2073      227731 :         an = hgmmodif(an0, PPOL);
    2074      227731 :         den = gmul(t0k, RgV_dotproduct(vecK, an));
    2075      227731 :         eno = gdiv(RgV_dotproduct(vecKi, an), den);
    2076      227731 :         if (gexpo(gsubgs(gabs(eno, prec), 1)) < -bitprec / 2)
    2077             :         {
    2078          98 :           if (pole)
    2079          77 :             for (i = 1; i < lB; i++)
    2080             :             {
    2081          42 :               GEN t = deg1pol_shallow(negi(powuu(BAD[i], k2)), gen_1, 0);
    2082          42 :               gmael(PPOL, i, 2) = gmul(gmael(PPOL, i, 2), t);
    2083             :             }
    2084         231 :           for (i = 1; i < lB; i++) gmael(PPOL, i, 2) = ginv(gmael(PPOL, i, 2));
    2085          98 :           eno = grndtoi(eno, NULL);
    2086          98 :           if (typ(eno) != t_INT) pari_err_BUG("lfunhgmwild");
    2087          98 :           return mkvec3(N, mkvec2(t, PPOL), eno);
    2088             :         }
    2089             :       }
    2090         756 :       set_avma(av0);
    2091             :     }
    2092             :   }
    2093             : }
    2094             : 
    2095             : static GEN
    2096          91 : BAD2small(GEN BAD)
    2097             : {
    2098          91 :   long i, l = lg(BAD);
    2099          91 :   GEN v = cgetg(l, t_VECSMALL);
    2100         210 :   for (i = 1; i < l; i++) v[i] = itou( gmael(BAD,i,1) );
    2101          91 :   return v;
    2102             : }
    2103             : 
    2104             : /* moments of motive */
    2105             : static GEN
    2106          91 : hgmmoments(GEN H, GEN t, GEN M, long nb)
    2107             : {
    2108          91 :   pari_sp av = avma, av2;
    2109          91 :   GEN P, v, S, L = hgmlfuninfty(H, t), BAD = BAD2small(gel(L, 3));
    2110          91 :   GEN k2 = gmul2n(gsubgs(gel(L,2), 1), -1);
    2111          91 :   long ct = 0, i, j, lP, lm, tm = typ(M);
    2112             : 
    2113          91 :   if (!nb) nb = 1000;
    2114          91 :   P = primes_zv(nb); lP = lg(P);
    2115          91 :   v = hgmcoefs(H, t, P[lP - 1]);
    2116          91 :   if (tm != t_VEC) M = gtovec(M);
    2117          91 :   av2 = avma; lm = lg(M);
    2118          91 :   S = const_vec(lm - 1, real_0(DEFAULTPREC));
    2119        3731 :   for (i = 1; i < lP; i++)
    2120             :   {
    2121        3640 :     long p = P[i];
    2122        3640 :     if (!Q_lval(t, p) && !zv_search(BAD, p))
    2123             :     {
    2124        3500 :       GEN T = gdiv(gel(v, p), gpow(utoipos(p), k2, DEFAULTPREC));
    2125        3500 :       ct++;
    2126        7000 :       for (j = 1; j < lm; j++)
    2127        3500 :         gel(S,j) = gadd(gel(S,j), gpow(T, gel(M,j), DEFAULTPREC));
    2128             :     }
    2129        3640 :     if ((i & 0xf) == 0) S = gerepilecopy(av2, S);
    2130             :   }
    2131          91 :   if (tm != t_VEC && tm != t_COL && tm != t_VECSMALL) S = gel(S, 1);
    2132          91 :   return gerepileupto(av, gdivgu(S, ct));
    2133             : }
    2134             : 
    2135             : /* Heuristic guess: is there a pole ? */
    2136             : static long
    2137          91 : lfunhgmispole(GEN H, GEN t, long nb)
    2138             : {
    2139          91 :   pari_sp av = avma;
    2140             :   GEN M;
    2141          91 :   if (!nb) nb = 40;
    2142          91 :   M = hgmmoments(H, t, gen_1, nb);
    2143          91 :   return gc_bool(av, gexpo(M) > -2);
    2144             : }
    2145             : 
    2146             : static GEN
    2147         112 : tag(GEN x, long t) { return mkvec2(mkvecsmall(t), x); }
    2148             : 
    2149             : static GEN
    2150         112 : lfunhgm_i(GEN hgm, GEN t, GEN hint, long bitprec)
    2151             : {
    2152         112 :   GEN L, vr, v = hgmlfuninfty(hgm, t), vga = zv_to_ZV(gel(v,1)), k = gel(v,2);
    2153         112 :   GEN BAD = gel(v,3), COND = gel(v, 4);
    2154         112 :   long pole = mpodd(k) && lfunhgmispole(hgm, t, 0), lB = lg(BAD);
    2155             : 
    2156         112 :   L = mkvecn(pole? 7: 6, tag(mkvec2(hgm,t), t_LFUN_HGM), gen_0, vga, k, COND,
    2157             :              gen_0, mkvec(mkvec2(shifti(addiu(k,1),-1), gen_0)));
    2158         112 :   if (pole && ldata_get_degree(L) == 1)
    2159             :   { /* motive = zeta */
    2160             :     long i;
    2161           0 :     gel(L,5) = gel(L,6) = gel(L,7) = gen_1;
    2162           0 :     if (lB != 1)
    2163             :     {
    2164           0 :       GEN R = mkrfrac(gen_1, deg1pol_shallow(gen_m1,gen_1,0));
    2165           0 :       gmael3(L, 1, 2, 2) = mkvec2(t, BAD);
    2166           0 :       for (i = 1; i < lB; i++) gmael(BAD, i, 2) = R;
    2167             :     }
    2168           0 :     return L;
    2169             :   }
    2170         112 :   if (lB == 1)
    2171             :   {
    2172          14 :     vr = lfunrootres(L, bitprec);
    2173          14 :     gel(L, 6) = gel(vr, 3);
    2174          14 :     if (pole) gel(L, 7) = gel(vr, 1);
    2175          14 :     return L;
    2176             :   }
    2177          98 :   v = lfunhgmwild(L, hgm, t, BAD, pole, hint, bitprec);
    2178          98 :   gel(L, 5) = gel(v, 1); /* N */
    2179          98 :   gmael3(L, 1, 2, 2) = gel(v, 2); /* [t, PPOL] */
    2180          98 :   gel(L, 6) = gel(v, 3); /* w */
    2181          98 :   if (pole)
    2182             :   {
    2183          35 :     vr = lfunrootres(L, bitprec);
    2184          35 :     gel(L, 7) = gel(vr, 1);
    2185             :   }
    2186          98 :   return L;
    2187             : }
    2188             : GEN
    2189         119 : lfunhgm(GEN hgm, GEN t, GEN hint, long bit)
    2190             : {
    2191         119 :   pari_sp av = avma;
    2192         119 :   if (!checkhgm(hgm)) pari_err_TYPE("lfunhgm", hgm);
    2193         112 :   return gerepilecopy(av, lfunhgm_i(hgm, t, hint, bit));
    2194             : }
    2195             : 
    2196             : GEN
    2197        3997 : dirhgm_worker(GEN P, ulong X, GEN hgm, GEN t)
    2198             : {
    2199        3997 :   pari_sp av = avma;
    2200        3997 :   long i, l = lg(P);
    2201        3997 :   GEN W = cgetg(l, t_VEC);
    2202       40822 :   for(i = 1; i < l; i++)
    2203             :   {
    2204       36825 :     ulong p = uel(P,i);
    2205       36825 :     long e, d = ulogint(X, p) + 1; /* minimal d such that p^d > X */
    2206       36825 :     gel(W,i) = RgXn_inv(hgmeulerfactorlimit(hgm, t, p, d-1, 0, &e), d);
    2207             :   }
    2208        3997 :   return gerepilecopy(av, mkvec2(P,W));
    2209             : }
    2210             : 
    2211             : GEN
    2212         399 : hgmcoefs(GEN hgm, GEN t, long n)
    2213             : {
    2214         399 :   GEN worker, bad = NULL;
    2215         399 :   if (!checkhgm(hgm)) pari_err_TYPE("hgmcoefs", hgm);
    2216         399 :   if (typ(t) == t_VEC && lg(t) == 3) { bad = gel(t,2); t = gel(t,1); }
    2217         399 :   if (!is_rational_t(typ(t))) pari_err_TYPE("hgmcoefs",t);
    2218         392 :   worker = snm_closure(is_entry("_dirhgm_worker"), mkvec2(hgm, t));
    2219         392 :   return pardireuler(worker, gen_2, stoi(n), NULL, bad);
    2220             : }

Generated by: LCOV version 1.16