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.14.0 lcov report (development 27775-aca467eab2) Lines: 1214 1253 96.9 %
Date: 2022-07-03 07:33:15 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      228222 : pari_init_hgm(void) { HGMDATA2 = HGMDATA3 = NULL; }
      25             : 
      26             : void
      27      226125 : pari_close_hgm(void)
      28      226125 : { 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      979398 : RgV_mirror(GEN x) { pari_APPLY_same(gsubsg(1, gel(x,i))); }
      46             : 
      47             : static GEN
      48      104139 : hodge(GEN val, GEN vbe, long *pTT)
      49             : {
      50      104139 :   long r = lg(val) - 1, T, s, k;
      51      104139 :   GEN H, w = indexsort( shallowconcat(val, RgV_mirror(vbe)) );
      52     1854657 :   for (k = 1, s = T = 0; k <= 2*r; k++)
      53     1750518 :     if (w[k] > r) { if (--s < T) T = s; } else s++;
      54      104139 :   H = zero_zv(r + 1 - T);
      55     1854657 :   for (k = 1, s = 0; k <= 2*r; k++) /* (x^-T * Hodge) += x^s, s - T >= 0 */
      56     1750518 :     if (w[k] > r) s--; else { H[s + 1 - T]++; s++; }
      57      104139 :   *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        1400 : al2cyE(GEN v)
      64             : {
      65             :   GEN dv, E, F, w;
      66        1400 :   long N, i, j, l = lg(v);
      67        1400 :   if (l == 1) return cgetg(1, t_VECSMALL);
      68        1232 :   w = Q_remove_denom(v, &dv);
      69        1232 :   if (!dv) return mkvecsmall(l-1);
      70         987 :   N = itou(dv); w = ZV_to_Flv(w, N); vecsmall_sort(w);
      71         987 :   E = zero_zv(N); /* to record polcyclo(i)^E[i] */
      72         987 :   F = cgetg(l, t_VECSMALL); /* to check input */
      73        4879 :   for (i = j = 1; i < l; i++)
      74             :   {
      75        3892 :     long k, m = w[i];
      76        3892 :     if (!m) { E[1]++; F[j++] = 0; }
      77        3521 :     else if (N % m == 0)
      78             :     {
      79        1890 :       long d = N/m, km;
      80        1890 :       E[d]++;
      81        8946 :       for (k = 1, km = m; k <= d; k++, km += m)
      82        7056 :         if (ugcd(k,d) == 1) F[j++] = km;
      83             :     }
      84             :   }
      85         987 :   setlg(F,j); vecsmall_sort(F);
      86         987 :   return gequal(w,F)? E: NULL;
      87             : }
      88             : 
      89             : static int
      90         700 : cyE_intersect(GEN a, GEN b)
      91             : {
      92         700 :   long i, l = minss(lg(a), lg(b));
      93        2653 :   for (i = 1; i < l; i++) if (a[i] && b[i]) return 1;
      94         693 :   return 0;
      95             : }
      96             : static GEN
      97         700 : get_CYCLOE(GEN val, GEN vbe)
      98             : {
      99         700 :   GEN Ea = al2cyE(val), Eb = al2cyE(vbe);
     100         700 :   if (!Ea || !Eb || cyE_intersect(Ea, Eb))
     101           7 :     pari_err_TYPE("hgminit [not a Q motive]", mkvec2(val,vbe));
     102         693 :   return mkvec2(Ea, Eb);
     103             : }
     104             : 
     105             : /* R[n/d] += mu(d) * r, for all d | n */
     106             : static void
     107       10472 : moebiusadd(GEN R, long n, long r)
     108             : {
     109       10472 :   if (r) {
     110        1764 :     GEN D = divisorsu_moebius(gel(factoru(n), 1));
     111        1764 :     long j, l = lg(D);
     112        1764 :     R[n] += r; /* d = 1 */
     113        3647 :     for (j = 2; j < l; j++) { long d = D[j]; R[n / labs(d)] += d < 0? -r: r; }
     114             :   }
     115       10472 : }
     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         693 : get_VPOLGA(GEN E)
     119             : {
     120         693 :   GEN Ea = gel(E,1), Eb = gel(E,2), R;
     121         693 :   long i, la = lg(Ea), lb = lg(Eb), n = maxuu(la, lb) - 1;
     122             :   pari_sp av;
     123         693 :   R = zero_zv(n); av = avma;
     124        7728 :   for (i = 1; i < la; i++) moebiusadd(R, i, Ea[i]);
     125        4130 :   for (i = 1; i < lb; i++) moebiusadd(R, i, -Eb[i]);
     126        5726 :   for (i = n; i; i--) if (R[i]) break;
     127         693 :   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        1204 : QV_normalize(GEN a)
     159             : {
     160        1204 :   long l = lg(a), i;
     161        1204 :   GEN v = cgetg(l, t_VEC);
     162        5908 :   for (i = 1; i < l; i++)
     163             :   {
     164        4718 :     GEN c = gel(a, i);
     165        4718 :     if (!is_rational_t(typ(c)))
     166          14 :       pari_err_TYPE("hgminit [not rational params]",c);
     167        4704 :     gel(v, i) = gfrac(c);
     168             :   }
     169        1190 :   return sort(v);
     170             : }
     171             : 
     172             : static GEN
     173        1344 : mangoldtexp(long n)
     174             : {
     175        1344 :   GEN D = divisorsu_moebius(gel(factoru(n),1)), P = gen_1;
     176        1344 :   long i, l = lg(D);
     177        4081 :   for (i = 1; i < l; i++)
     178             :   { /* prod (n/d) ^ (n/d mu(d)) */
     179        2737 :     long mud = D[i], nd = n / mud; /* mud = mu(d) * d */
     180        2737 :     P = gmul(P, powis(utoi(labs(nd)), nd));
     181             :   }
     182        1344 :   return P;
     183             : }
     184             : static GEN
     185         924 : E2exp(GEN E)
     186             : {
     187         924 :   long n, r, l = lg(E);
     188         924 :   GEN P = gen_1;
     189        9786 :   for (n = 1; n < l; n++)
     190        8862 :     if ((r = E[n])) P = gmul(P, gpowgs(mangoldtexp(n), r));
     191         924 :   return P;
     192             : }
     193             : static long
     194         959 : 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         504 : initab(GEN a, GEN b)
     200             : {
     201             :   GEN VALNUM, VALDEN, VBENUM, VBEDEN, CYCLOE, SIGNPAR;
     202             :   GEN BAD, HODGE, VPOLGA, res;
     203         504 :   long j, WT, TT, OFFMPOL, l = lg(a), DEG = l-1, SWAP = 0;
     204             : 
     205         504 :   if (lg(b) != l) pari_err_TYPE("hgminit [#al != #be]", mkvec2(a,b));
     206         504 :   if (l == 1) pari_err_TYPE("hgminit [#al = 0]", mkvec2(a,b));
     207         490 :   a = QV_normalize(a);
     208         483 :   b = QV_normalize(b);
     209         476 :   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         469 :   VALNUM = cgetg(l, t_VECSMALL);
     216         469 :   VALDEN = cgetg(l, t_VECSMALL);
     217         469 :   VBENUM = cgetg(l, t_VECSMALL);
     218         469 :   VBEDEN = cgetg(l, t_VECSMALL);
     219        2443 :   for (j = 1; j < l; j++)
     220             :   {
     221        1974 :     Qtoss(gel(a,j), &VALNUM[j], &VALDEN[j]);
     222        1974 :     Qtoss(gel(b,j), &VBENUM[j], &VBEDEN[j]);
     223             :   }
     224         469 :   BAD = negi(lcmii(ZV_lcm(zv_to_ZV(VALDEN)), ZV_lcm(zv_to_ZV(VBEDEN))));
     225         469 :   HODGE = hodge(a, b, &TT);
     226         469 :   WT = degpol(HODGE);
     227         469 :   CYCLOE = get_CYCLOE(a, b);
     228         462 :   VPOLGA = get_VPOLGA(CYCLOE);
     229         462 :   SIGNPAR = odd(WT)? gen_1: discprod(gel(CYCLOE, odd(DEG)? 2: 1));
     230         462 :   OFFMPOL = (get_b1(CYCLOE) - zv_sum(VPOLGA)) >> 1;
     231         462 :   res = cgetg(13, t_VEC);
     232         462 :   gel(res, 1) = a;
     233         462 :   gel(res, 2) = b;
     234         462 :   gel(res, 3) = CYCLOE;
     235         462 :   gel(res, 4) = VBENUM;
     236         462 :   gel(res, 5) = VBEDEN;
     237         462 :   gel(res, 6) = gdiv(E2exp(gel(CYCLOE,1)), E2exp(gel(CYCLOE,2))); /*MVALUE*/
     238         462 :   gel(res, 7) = VPOLGA;
     239         462 :   gel(res, 8) = BAD;
     240         462 :   gel(res, 9) = HODGE;
     241         462 :   gel(res, 10) = get_u(a, b, CYCLOE, VPOLGA, DEG, WT);
     242         462 :   gel(res, 11) = SIGNPAR;
     243         462 :   gel(res, 12) = mkvecsmall3(OFFMPOL, TT, SWAP);
     244         462 :   return res;
     245             : }
     246             : static int
     247        1253 : checkhgm(GEN v)
     248             : {
     249        1246 :   return typ(v) == t_VEC && lg(v) == 13
     250        2499 :       && typ(gel(v,12)) == t_VECSMALL && lg(gel(v,12)) == 4;
     251             : }
     252             : 
     253             : /* accessors */
     254      141726 : 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        1582 : static GEN hgm_get_CYCLOE(GEN v) { return gel(v, 3); }
     257       37593 : static GEN hgm_get_VBENUM(GEN v) { return gel(v, 4); }
     258       37593 : static GEN hgm_get_VBEDEN(GEN v) { return gel(v, 5); }
     259         623 : static GEN hgm_get_MVALUE(GEN v) { return gel(v, 6); }
     260      687866 : static GEN hgm_get_VPOLGA(GEN v) { return gel(v, 7); }
     261       37827 : static GEN hgm_get_BAD(GEN v) { return gel(v, 8); }
     262      109707 : 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       28624 : static GEN hgm_get_SIGNPAR(GEN v) { return gel(v, 11); }
     266      141691 : static long hgm_get_DEG(GEN v) { return lg(hgm_get_VAL(v))-1; }
     267      652486 : static long hgm_get_OFFMPOL(GEN v) { return gel(v, 12)[1]; }
     268       76957 : static long hgm_get_TT(GEN v) { return gel(v, 12)[2]; }
     269         672 : static long hgm_get_SWAP(GEN v) { return gel(v, 12)[3]; }
     270      109343 : 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       59224 :   for (i = 2; i < p; i++) gel(v, i+1) = gdivgu(gel(v, i), i);
     285        2155 :   gel(w, 1) = shallowcopy(v);
     286      260485 :   for (k = s = 1, t = p; k <= L; k++)
     287             :   {
     288      258320 :     gel(v, s) = gdivgu(gadd(gel(v, t), gel(v, s)), k + p - 1);
     289      258330 :     if (++s > p) { s = 1; gel(w, k/p + 1) = shallowcopy(v); }
     290      258330 :     if (++t > p) t = 1;
     291             :   }
     292        2165 :   return w;
     293             : }
     294             : 
     295             : /* [binom(x + k - 1, k) k! p^k, k = 1..N], x a t_PADIC */
     296             : static GEN
     297      116943 : binomfact(GEN x, long N)
     298             : {
     299      116943 :   GEN z, v = cgetg(N+1, t_VEC);
     300             :   long i;
     301      116943 :   gel(v,1) = z = leafcopy(x);
     302      116944 :   setvalp(z, valp(z)+1);
     303      625552 :   for (i = 1; i < N; i++)
     304             :   {
     305      508629 :     gel(v, i+1) = z = gmul(z, gaddgs(x, i));
     306      508557 :     setvalp(z, valp(z) + 1);
     307             :   }
     308      116923 :   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      116943 : sumbinom(GEN v, GEN x, long m0, long N, ulong pD)
     313             : {
     314      116943 :   pari_sp av = avma;
     315      116943 :   GEN C = binomfact(x, N), S = gmael(v, 1, m0+1); /* k = 0 */
     316             :   long k;
     317      742408 :   for (k = 1; k <= N; k++) S = gadd(S, gmul(gel(C,k), gmael(v, k+1, m0+1)));
     318      116923 :   return gc_ulong(av, Rg_to_Fl(S, pD));
     319             : }
     320             : /* x * y mod p^D, y a p-unit */
     321             : static GEN
     322      203890 : umultop(ulong x, ulong y, GEN gp, GEN gpD, long D)
     323             : {
     324      203890 :   ulong pD = gpD[2];
     325             :   long v;
     326             :   GEN z;
     327      203890 :   if (!x) return zeropadic_shallow(gp, D);
     328      203890 :   v = u_lvalrem(x, gp[2], &x);
     329      203890 :   if (x >= pD) x %= pD;
     330      203890 :   z = cgetg(5, t_PADIC);
     331      203889 :   z[1] = evalprecp(D) | evalvalp(v);
     332      203888 :   gel(z,2) = gp;
     333      203888 :   gel(z,3) = gpD;
     334      203888 :   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      104477 :   for (m = 1; m <= q; m += 2, set_avma(av2))
     368             :   {
     369      102456 :     if (f > 1) { m0 = m%p; m1 = m / p; M = m1 + m0 * pf1; } else M = m0 = m;
     370      102456 :     x = umultop(M, iQ, gp, gpD, D);
     371      102453 :     vres[m+1] = S = sumbinom(vall, x, m0, N, pD);
     372      102450 :     if (!odd(m0)) S = pD - S;
     373      102450 :     vres[pf - m] = Fl_inv(S, pD);
     374             :   }
     375       16458 :   for (m = vpowv; m <= q; m += 2*vpowv, set_avma(av2))
     376             :   {
     377       14437 :     if (f > 1) { m0 = m%p; m1 = m / p; M = m1 + m0 * pf1; } else M = m0 = m;
     378       14437 :     x = umultop(M, iQ, gp, gpD, D);
     379       14438 :     vres[m+1] = S = sumbinom(vall, x, m0, N, pD);
     380       14437 :     if (!odd(m0)) S = pD - S;
     381       14437 :     vres[pf - m] = Fl_inv(S, pD);
     382             :   }
     383        2021 :   l2 = gmulgu(Qp_log(cvtop(gen_2, gp, D)), p - 1);
     384        2020 :   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       93796 :     for (m = vpo; m <= q; m += 2*vpo, set_avma(av2))
     390             :     {
     391       86946 :       ulong ms2 = (m >> 1) + 1;
     392       86946 :       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       86946 :       m0 = p - m%p; m1 = m/p + 1; M = pM + m1 - m0 * pf1;
     396       86946 :       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       90433 :   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       33935 : teichmodp2(GEN VPOLGA, ulong p, ulong p2)
     475             : {
     476             :   ulong pN;
     477             :   long c, l, N;
     478       33935 :   GEN W = cgetg_copy(VPOLGA,&l); /* W[N] = teich(N)^N */
     479      142852 :   for (N = 1, pN = p; N < l; N++, pN += p)
     480      108923 :     if ((c = VPOLGA[N]))
     481             :     {
     482       82247 :       long N0; (void)z_lvalrem(N, p, &N0);
     483       82246 :       if (c < 0) N0 = Fl_inv(N0 % p, p);
     484       82245 :       W[N] = Fl_powu(N0, pN, p2);
     485             :     }
     486       33929 :   return W;
     487             : }
     488             : /* GF[i] = i! mod p^2 */
     489             : static GEN
     490       33930 : get_GF(ulong p, ulong p2)
     491             : {
     492       33930 :   GEN F = cgetg(p, t_VECSMALL);
     493             :   ulong i;
     494    33022197 :   F[1] = 1; F[2] = 2; for (i = 3; i < p; i++) F[i] = Fl_mul(F[i-1], i, p2);
     495       33068 :   return F;
     496             : }
     497             : 
     498             : /* p odd, return v s.t. v[i] = 1/i mod p, i < p/2 */
     499             : static GEN
     500       33932 : Flv_inv_p2(ulong p)
     501             : {
     502       33932 :   long i, g = pgener_Fl(p), ph = (p >> 1);
     503       33934 :   GEN v = cgetg(ph + 1, t_VECSMALL);
     504       33933 :   pari_sp av = avma;
     505       33933 :   GEN w = cgetg(ph, t_VECSMALL);
     506       34443 :   w[1] = g;
     507    16641230 :   for (i = 2; i < ph; i++) w[i] = Fl_mul(w[i-1], g, p); /* w[i]=g^i */
     508    17367665 :   for (i = 1; i < ph; i++)
     509             :   {
     510    17333742 :     long x = w[i], y = w[ph - i]; /* g^(-i) = - g^(ph - i) */
     511    17333742 :     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       33926 : get_GH(GEN F, ulong p)
     519             : {
     520       33926 :   long i, g, ph = p >> 1;
     521       33926 :   GEN H = cgetg(p, t_VECSMALL), v = Flv_inv_p2(p);
     522       33936 :   H[p-1] = (F[p-1] + 1) / p; g = p - H[p-1];
     523    16184200 :   for (i = 1; i <= ph; i++)
     524             :   {
     525    16150282 :     long c = g = Fl_add(g, v[i], p);
     526    16144175 :     H[i] = c = Fl_mul(c, i, p);
     527    16152033 :     H[p-1 - i] = Fl_neg(Fl_add(c, g, p), p);
     528             :   }
     529       33918 :   return H;
     530             : }
     531             : /* p odd, GPB[m+1] = Gamma_p(m/(p-1)) mod p^2 */
     532             : static GEN
     533       33929 : doprecompmodp2(ulong p, ulong p2)
     534             : {
     535       33929 :   GEN F = get_GF(p, p2), H = get_GH(F, p), GPV = cgetg(p, t_VECSMALL);
     536             :   ulong r;
     537    32316993 :   for (r = 1; r < p; r++)
     538             :   {
     539    32283081 :     long t = Fl_mul(F[r], 1 + p * H[r], p2);
     540    32273841 :     GPV[p - r] = odd(r)? t: p2 - t;
     541             :   }
     542       33912 :   return GPV;
     543             : }
     544             : 
     545             : /***************************************************************/
     546             : /*        PART III: Motive initializations depending on p      */
     547             : /***************************************************************/
     548             : static GEN
     549       37593 : get_L1(GEN hgm, long PFM1, long f)
     550             : {
     551       37593 :   GEN VBEDEN = hgm_get_VBEDEN(hgm);
     552       37593 :   GEN VBENUM = hgm_get_VBENUM(hgm), v;
     553       37593 :   long DEG = hgm_get_DEG(hgm), j;
     554       37593 :   v = const_vecsmall(PFM1, f * hgm_get_TT(hgm));
     555      164019 :   for (j = 1; j <= DEG; j++)
     556      126423 :     if (PFM1 % VBEDEN[j] == 0)
     557             :     {
     558      118301 :       long t = ((-(PFM1 / VBEDEN[j]) * VBENUM[j]) % PFM1) + 1;
     559      118301 :       if (t <= 0) t += PFM1;
     560      118301 :       v[t] -= f;
     561             :     }
     562       37596 :   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       33963 : get_L0(GEN hgm, ulong PM1)
     569             : {
     570       33963 :   GEN perm, vt, vc, VL0, VPOLGA = hgm_get_VPOLGA(hgm);
     571       33963 :   long w, S, c, d, j, N, m, l = lg(VPOLGA), D = (l * (l-1)) >> 1;
     572       33963 :   vt = cgetg(D+1, t_VECSMALL);
     573       33966 :   vc = cgetg(D+1, t_VECSMALL);
     574      109086 :   for (N = 2, c = 1; N < l; N++)
     575       75123 :     if (VPOLGA[N])
     576             :     {
     577             :       ulong Q;
     578             :       long q;
     579      195462 :       for (q = 1, Q = 0; q <= N; q++, Q += PM1, c++)
     580             :       {
     581      147097 :         vt[c] = ceildivuu(Q, N);
     582      147092 :         vc[c] = VPOLGA[N];
     583             :       }
     584             :     }
     585       33963 :   setlg(vt,c); setlg(vc,c); perm = vecsmall_indexsort(vt);
     586       33964 :   vt = vecsmallpermute(vt, perm);
     587       33963 :   vc = vecsmallpermute(vc, perm);
     588       33963 :   w = vt[1];
     589      147138 :   for (j = 2, d = 1; j < c; j++)
     590      113175 :     if (vt[j] == w) vc[d] += vc[j];
     591             :     else
     592             :     {
     593       79288 :       d++;
     594       79288 :       vt[d] = w = vt[j];
     595       79288 :       vc[d] = vc[j];
     596             :     }
     597       33963 :   d++; vt[d] = PM1; vc[d] = 0; /* sentinels */
     598       33963 :   VL0 = cgetg(PM1+1, t_VECSMALL); S = hgm_get_OFFMPOL(hgm);
     599      147207 :   for (j = 1; j < d; j++)
     600             :   {
     601    34875919 :     for (m = vt[j]; m < vt[j+1]; m++) VL0[m+1] = S;
     602      113240 :     S -= vc[j+1];
     603             :   }
     604       33967 :   return VL0;
     605             : }
     606             : /* 2) direct computation (single m), good when p << N log N */
     607             : static long
     608      614864 : L0(GEN hgm, ulong p, long f, long PFM1, long m)
     609             : {
     610      614864 :   GEN VPOLGA = hgm_get_VPOLGA(hgm);
     611      614864 :   long e, S = hgm_get_OFFMPOL(hgm) * f, l = lg(VPOLGA);
     612             :   ulong pem;
     613     1858432 :   for (e = 0, pem = m; e < f; e++, pem *= p)
     614             :   {
     615     1243569 :     long N, s = 0;
     616             :     ulong Npem;
     617     9691508 :     for (N = 1, Npem = pem; N < l; N++, Npem += pem)
     618     8447939 :       if (VPOLGA[N]) s += VPOLGA[N] * (long)(Npem / PFM1);
     619     1243569 :     S -= s;
     620             :   }
     621      614863 :   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      984441 : gapnpow(GEN T, long p, long f, long PFM1, long N, ulong Nmpe)
     656             : {
     657             :   GEN Vr;
     658             :   long Nm, i, S;
     659             : 
     660      984441 :   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      235072 : hgmC(GEN VPOLGA, GEN GPV, GEN TEICH, ulong p, long f, ulong PFM1, ulong m, long dfp)
     673             : {
     674      235072 :   GEN r = gen_1;
     675      235072 :   long c, e, N, l = lg(VPOLGA);
     676             :   ulong Nmpe, mpe;
     677      700153 :   for (e = 0, mpe = m; e < f; e++, mpe = Fl_mul(mpe, p, PFM1))
     678     3819361 :     for (N = 1, Nmpe = mpe; N < l; N++, Nmpe = Fl_add(Nmpe, mpe, PFM1))
     679     3354278 :       if ((c = VPOLGA[N]))
     680             :       { /* Gamma_p(frac(N*m*p^e/(p^f-1)))  */
     681     1449530 :         GEN z = utoi(GPV[Nmpe + 1]);
     682     1449535 :         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     1449512 :         r = gmul(r, gpowgs(z, c));
     685             :       }
     686      235054 :   if (typ(r) != t_PADIC) r = cvtop(r, utoipos(p), dfp);
     687      235053 :   return r;
     688             : }
     689             : 
     690             : /* Same modulo p^2, Wm[N] = teich(N^(sign(VPOLGA[N]) * m)) */
     691             : static ulong
     692    16405384 : hgmCmodp2(GEN VPOLGA, GEN Wm, GEN GPV, ulong PFM1, ulong m, ulong p2)
     693             : {
     694    16405384 :   long l = lg(VPOLGA), c, N;
     695    16405384 :   ulong res = 1, Nm;
     696    64678005 :   for (N = 1, Nm = m; N < l; N++, Nm = Fl_add(Nm, m, PFM1))
     697    48335900 :     if ((c = VPOLGA[N]))
     698             :     {
     699             :       ulong z;
     700    36981506 :       if (c > 0) z = GPV[Nm + 1];
     701             :       else
     702             :       {
     703    19008045 :         c = -c;
     704    19008045 :         if (!Nm) z = 1;
     705             :         else
     706             :         {
     707    18988652 :           z = GPV[PFM1 + 1 - Nm];
     708    18988652 :           if (!odd(Nm)) z = p2 - z; /* GPV[Nm+1]^(-1) by reflection formula */
     709             :         }
     710             :       }
     711    37672769 :       if (N > 1) z = Fl_mul(z, Wm[N], p2);
     712    37619497 :       z = Fl_powu(z, c, p2);
     713             :       /* z = (GPV[Nm+1] * teich(N^m))^VPOLGA[N] */
     714    37032607 :       res = res == 1? z: Fl_mul(res, z, p2);
     715             :     }
     716    16266251 :   return res;
     717             : }
     718             : 
     719             : /***************************************************************/
     720             : /*        PART IV: Motive functions depending on p, t          */
     721             : /***************************************************************/
     722             : static long
     723       37833 : get_dfp(GEN hgm, long p, long f)
     724             : {
     725       37833 :   long DEG = hgm_get_DEG(hgm), WT = hgm_get_WT(hgm);
     726       37833 :   long dfp = ceil(log((double)2*DEG) / log((double)p) + f * (WT * 0.5 ) + 1e-5);
     727       37833 :   if (p == 2) dfp++; /* FIXME: why ? */
     728       37833 :   return dfp;
     729             : }
     730             : /* V being a vecsmall, compute all p^TT*hgmC(m)/hgmC(0) for indices in V */
     731             : static GEN
     732        3667 : hgmCall(GEN hgm, ulong p, long f, long dfp, GEN V)
     733             : {
     734        3667 :   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      618600 :   for (; i < lV; i++)
     756             :   {
     757      614940 :     long m = V? V[i]: i-1;
     758      614940 :     long L = VL0? VL0[m+1]: L0(hgm, p, f, PFM1, m);
     759      614940 :     long e = L + VL1[m+1];
     760      614940 :     if (!V && e >= dfp) c = gen_0;
     761             :     else
     762             :     { /* FIXME: dfp is fishy in Jordantame, don't trust if V = NULL */
     763      235071 :       pari_sp av = avma;
     764      235071 :       c = hgmC(VPOLGA, GPV, TEICH, p, f, PFM1, m, dfp);
     765      235053 :       setvalp(c, e); if (odd(L ^ l0)) c = gneg(c);
     766      235050 :       c = gerepileupto(av, padic_to_Q(c));
     767             :     }
     768      614939 :     gel(v, V? i: i+1) = c;
     769             :   }
     770        3660 :   return v;
     771             : }
     772             : /* Same mod p^2, f = 1 */
     773             : static GEN
     774       33936 : hgmCallmodp2(GEN hgm, ulong p)
     775             : {
     776       33936 :   GEN C, GPV, VL0, VL1, W1, Wm, VPOLGA = hgm_get_VPOLGA(hgm);
     777       33935 :   long l = lg(VPOLGA), TT = hgm_get_TT(hgm);
     778       33935 :   ulong m, PFM1 = p - 1, p2 = p * p;
     779             : 
     780       33935 :   if (p & HIGHMASK) pari_err_OVERFLOW("hgmCallmodp2");
     781       33928 :   VL0 = get_L0(hgm, PFM1);
     782       33932 :   VL1 = get_L1(hgm, PFM1, 1);
     783       33935 :   W1 = teichmodp2(VPOLGA, p, p2);
     784       33930 :   Wm = shallowcopy(W1);
     785       33929 :   GPV = doprecompmodp2(p, p2);
     786       33913 :   C = cgetg(PFM1+2, t_VECSMALL); C[1] = evalvarn(0);
     787       33923 :   C[2] = TT > 1? 0: (TT == 1? p : 1); /* m = 0 */
     788    31789548 :   for (m = 1; m < PFM1; m++)
     789             :   {
     790    32414490 :     long e = VL0[m + 1] + VL1[m + 1], j;
     791             :     ulong c;
     792    32414490 :     if (e >= 2) c = 0;
     793             :     else
     794             :     {
     795    16325001 :       c = hgmCmodp2(VPOLGA, Wm, GPV, PFM1, m, p2);
     796    16261309 :       if (odd(VL0[m + 1] ^ VL0[1])) c = Fl_neg(c, p2);
     797    15675907 :       if (e == 1) c = (c % p) * p;
     798             :     }
     799    31765396 :     C[m + 2] = c;
     800    94353382 :     for (j = 2; j < l; j++)
     801    62597757 :       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        3850 : inv(ulong q, long n)
     809             : {
     810        3850 :   ulong z = q + 1;
     811             :   long i;
     812        7980 :   for (i = 1; i < n; i++) z = z * q + 1;
     813        3850 :   return z;
     814             : }
     815             : 
     816             : /* General H function: C(teich(f + O(p^dfp))) / (1 - p^f) */
     817             : static GEN
     818        3849 : hgmH(GEN C, long p, long f, long dfp, GEN t)
     819             : {
     820        3849 :   GEN q = powuu(p, dfp), z;
     821             :   long n;
     822        3849 :   (void)Q_lvalrem(t, p, &t);
     823        3849 :   z = Rg_to_Fp(t, q);
     824        3849 :   z = Zp_teichmuller(z, utoipos(p), dfp, q);
     825        3849 :   z = FpX_eval(C, z, q);
     826        3850 :   n = dfp / f; if (!(dfp % f)) n--;
     827        3850 :   z = Fp_mulu(z, inv(upowuu(p,f), n), q);
     828        3850 :   return Fp_center(z, q,  shifti(q,-1));
     829             : }
     830             : static GEN
     831       33922 : hgmHmodp2(GEN C, ulong p, GEN t)
     832             : {
     833       33922 :   ulong p2 = p * p, wt = Fl_powu(Rg_to_Fl(t, p2), p, p2);
     834       33934 :   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       37791 : hgmtrace(GEN hgm, long p, long f, GEN t, long c)
     842             : {
     843       37791 :   long dfp = get_dfp(hgm, p, f);
     844       37792 :   if (c == C_FAKE) return hgmU(hgm, p, f, t, dfp);
     845       37561 :   if (f == 1 && dfp <= 2) return hgmHmodp2(hgmCallmodp2(hgm, p), p, t);
     846        3625 :   return hgmH(hgmCall(hgm, p, f, dfp, NULL), p, f, dfp, t);
     847             : }
     848             : 
     849             : static GEN
     850         693 : F2v_factorback(GEN E)
     851             : {
     852         693 :   long i, l = lg(E);
     853         693 :   GEN c = gen_1;
     854        4186 :   for (i = 1; i < l; i++)
     855        3493 :     if (odd(E[i])) c = muliu(c, i);
     856         693 :   return c;
     857             : }
     858             : 
     859             : static long
     860       28414 : Q_krois(GEN T, long p) { return krouu(Rg_to_Fl(T, (p == 2)? 8: p), p); }
     861             : static long
     862       33926 : hgmsign(GEN hgm, long p, GEN t)
     863             : {
     864             :   GEN u;
     865       33926 :   if (odd(hgm_get_WT(hgm))) return 1;
     866       28414 :   t = ginv(t); u = gmul(gsubsg(1, t), hgm_get_SIGNPAR(hgm));
     867       28413 :   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         245 : hgmcharcond(GEN hgm, GEN t)
     872             : {
     873             :   GEN u;
     874         245 :   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         364 : upowers_u(ulong r, long n, ulong a)
     893             : {
     894         364 :   long i, l = n + 2, t = a;
     895         364 :   GEN v = cgetg(l, t_VECSMALL);
     896        2310 :   for (i = 1; i < l; i++) { v[i] = t; t *= r; }
     897         364 :   return v;
     898             : }
     899             : static GEN
     900       36871 : powers_u(ulong r, long n, ulong a)
     901             : {
     902       36871 :   long i, l = n + 2;
     903       36871 :   GEN v = cgetg(l, t_VEC), t = utoi(a);
     904      106574 :   for (i = 1; i < l; i++) { gel(v,i) = t; t = muliu(t, r); }
     905       36870 :   return v;
     906             : }
     907             : static GEN
     908       36866 : mkpowers(long p, long d, long WT)
     909             : {
     910             :   ulong q, r;
     911       36866 :   if (WT == 0) q = r = 1;
     912       35690 :   else if (!odd(d)) q = r = upowuu(p, WT);
     913       21886 :   else if (WT == 1) { q = 1; r = p; }
     914       21886 :   else { q = upowuu(p, WT >> 1); r = q*q; if (odd(WT)) r *= p; }
     915       36869 :   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       36868 : Efuneq(GEN E, long p, long d, long WT, long SIGN, long B)
     922             : {
     923       36868 :   long j = maxss(d - B, 0), k = d + 1 - j, nE = lg(E)-1, l = (d + 1) >> 1;
     924       36867 :   GEN T = cgetg(k + 1, t_VEC), vp = mkpowers(p, d, WT);
     925       39705 :   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       39494 :   for (; k >= nE; k--) gel(T, k) = gen_0;
     932      108595 :   for (; k > 0; k--) gel(T, k) = gel(E, k+1);
     933       36870 :   return RgV_to_RgX(T,0);
     934             : }
     935             : 
     936             : static long
     937       37624 : hgmclass(GEN hgm, long p, GEN t)
     938             : {
     939       37624 :   GEN BAD = hgm_get_BAD(hgm);
     940             :   long ap, bp;
     941       37624 :   if (!umodiu(BAD, p))
     942             :   {
     943         665 :     long v = Q_lval(t, p);
     944         665 :     if (v && v + Q_lval(hgm_get_MVALUE(hgm), p) == 0)
     945             :     {
     946         147 :       GEN N = hgmcharcond(hgm, t);
     947         147 :       if (umodiu(N, p)) return C_FAKE; /* a(chi) = 0 */
     948             :     }
     949         553 :     return C_BAD;
     950             :   }
     951       36960 :   if (typ(t) == t_INT)
     952             :   {
     953       23992 :     ap = umodiu(t, p); if (!ap) return C_TAME0;
     954       23965 :     bp = 1;
     955             :   }
     956             :   else
     957             :   {
     958       12968 :     ap = umodiu(gel(t,1), p); if (!ap) return C_TAME0;
     959       12962 :     bp = umodiu(gel(t,2), p); if (!bp) return C_TAME0;
     960             :   }
     961       36864 :   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       36947 : frobpoltrunc(GEN hgm, GEN t, long c, long p, long B, long *pF)
     968             : {
     969       36947 :   long DEG = hgm_get_DEG(hgm), WT = hgm_get_WT(hgm);
     970       36946 :   long D = isint1(t)? (odd(WT) ? DEG - 2 : DEG - 1): DEG;
     971       36945 :   long f, mi, m, q = upowuu(p, WT >> 1);
     972             :   GEN E, s;
     973             : 
     974       36946 :   mi = minss(B, (c == C_FAKE)? D: D >> 1);
     975       36945 :   m = (mi == D && c == C_TAME1 && !odd(DEG))? mi: mi+1;
     976       36945 :   s = cgetg(m+1, t_POL); s[1] = evalsigne(1)|evalvarn(0);
     977       74723 :   for (f = 1; f < m; f++) gel(s, f+1) = negi( hgmtrace(hgm, p, f, t, c) );
     978       36937 :   *pF = 0; s = RgXn_expint(s, m);
     979       36931 :   if (mi == D) return s;
     980       36868 :   if (c == C_TAME1)
     981             :   {
     982        2940 :     long SIGN = kroiu(hgm_get_U(hgm), p);
     983        2940 :     if (odd(DEG))
     984         203 :       E = Efuneq(s, p, DEG-1, WT, SIGN, B);
     985             :     else
     986             :     {
     987        2737 :       GEN T = deg1pol_shallow(stoi(- SIGN * q), gen_1, 0);
     988        2737 :       E = RgXn_mul(s, RgXn_inv(T, m), m);
     989        2737 :       E = Efuneq(E, p, DEG - 2, WT, 1, B);
     990        2737 :       if (!gequal1(t) || !odd(WT)) E = gmul(E, T);
     991             :     }
     992        2940 :     *pF = 1;
     993        2940 :     if (!odd(WT))
     994             :     {
     995             :       GEN T, u, t0;
     996         217 :       long v = Q_lvalrem(gsubgs(t, 1), p, &t0);
     997         217 :       if (!odd(v))
     998             :       {
     999          63 :         if (typ(t0) == t_FRAC) t0 = mulii(gel(t0,1), gel(t0,2));
    1000          63 :         u = coredisc(mulii(t0, hgm_get_U0(hgm)));
    1001          63 :         T = deg1pol_shallow(stoi(-kroiu(u, p) * q), gen_1, 0);
    1002          63 :         E = gmul(E, T); *pF = 0;
    1003             :       }
    1004             :     }
    1005             :   }
    1006             :   else
    1007             :   {
    1008       33928 :     long SIGN = hgmsign(hgm, p, t);
    1009       33928 :     E = Efuneq(s, p, DEG, WT, SIGN, B);
    1010             :   }
    1011       36869 :   return E;
    1012             : }
    1013             : 
    1014             : GEN
    1015          84 : hgmcoef(GEN hgm, GEN t, GEN n)
    1016             : {
    1017          84 :   pari_sp av = avma;
    1018          84 :   GEN T, P, E, F = check_arith_all(n, "hgmcoef");
    1019             :   long i, lP;
    1020             : 
    1021          84 :   if (!checkhgm(hgm)) pari_err_TYPE("hgmcoef", hgm);
    1022          84 :   if (!is_rational_t(typ(t))) pari_err_TYPE("hgmcoef",t);
    1023          77 :   if (hgm_get_SWAP(hgm)) t = ginv(t);
    1024          77 :   if (!F) { F = Z_factor(n); P = gel(F,1); }
    1025             :   else
    1026             :   {
    1027          21 :     P = gel(F,1);
    1028          21 :     if (lg(P) == 1 || signe(gel(P,1)) <= 0) return gen_0;
    1029          21 :     n = typ(n) == t_VEC? gel(n,1): factorback(F);
    1030             :   }
    1031          77 :   if (signe(n) <= 0) pari_err_DOMAIN("hgmcoef", "n", "<=", gen_0, n);
    1032          77 :   E = gel(F,2); lP = lg(P); T = gen_1;
    1033          77 :   T = gen_1;
    1034         161 :   for (i = 1; i < lP; i++)
    1035             :   {
    1036          84 :     long e, p = itos(gel(P, i)), f = itos(gel(E, i)), c = hgmclass(hgm, p, t);
    1037             :     GEN A;
    1038          84 :     if (c == C_BAD) pari_err_IMPL("hgmcoef for bad primes");
    1039          84 :     A = frobpoltrunc(hgm, t, c, p, f, &e);
    1040          84 :     T = gmul(T, RgX_coeff(RgXn_inv(A, f+1), f));
    1041             :   }
    1042          77 :   return gerepilecopy(av, T);
    1043             : }
    1044             : 
    1045             : static GEN
    1046          14 : count2list(GEN E)
    1047             : {
    1048          14 :   long i, j, r, l = lg(E);
    1049          14 :   GEN v = cgetg(zv_sum(E)+1, t_VECSMALL);
    1050          49 :   for (i = j = 1; i < l; i++)
    1051          56 :     if ((r = E[i])) while(r--) v[j++] = i;
    1052          14 :   setlg(v,j); return v;
    1053             : }
    1054             : static GEN hgminit_i(GEN val, GEN vbe);
    1055             : #if 0
    1056             : static long
    1057             : minval(GEN F, long p)
    1058             : {
    1059             :   long i, d = degpol(F), m = LONG_MAX;
    1060             :   for (i = 1; i <= d; i++) m = minss(m, Z_lval(gel(F,i+2), p) / i);
    1061             :   return m;
    1062             : }
    1063             : static GEN
    1064             : cycloE_filter(GEN A, long p)
    1065             : {
    1066             :   long i, j, l = lg(A);
    1067             :   GEN B = shallowcopy(A);
    1068             :   if (p < l)
    1069             :     for (i = p; i < l; i += p)
    1070             :       if (A[i])
    1071             :       {
    1072             :         (void)z_lvalrem(i / p, p, &j);
    1073             :         B[j] += A[i];
    1074             :         B[i] = 0;
    1075             :       }
    1076             :   return B;
    1077             : }
    1078             : /* doesn't work because A/B is not valid hypergeometric data */
    1079             : static GEN
    1080             : eulfacbadnew(GEN hgm, GEN t, long p, long *pe)
    1081             : {
    1082             :   GEN AB = hgm_get_CYCLOE(hgm), gp = utoipos(p);
    1083             :   GEN A = cycloE_filter(gel(AB,1), p);
    1084             :   GEN B = cycloE_filter(gel(AB,2), p);
    1085             :   GEN hgmp = hgminit_i(count2list(A), count2list(B));
    1086             :   GEN cF, E, F = hgmeulerfactor(hgmp, t, p, &E);
    1087             :   long s, v, d = degpol(F), w = hgm_get_WT(hgm);
    1088             :   *pe = 0;
    1089             :   if (!d) return pol_1(0);
    1090             :   cF = leading_coeff(F);
    1091             :   v = Z_lvalrem(cF, p, &cF);
    1092             :   if (!is_pm1(cF)) return pol_1(0);
    1093             :   s = minval(F, p);
    1094             :   if (s) { F = RgX_unscale(F, powis(gp, -s)); v -= s * d; }
    1095             :   if ((2 * v) % d || (!odd(w) && v % d)) return F;
    1096             :   return RgX_unscale(F, powis(gp, w / 2 - v / d));
    1097             : }
    1098             : #endif
    1099             : static GEN Jordantameexpo(GEN hgm, long v, GEN t0, long p, long *pe);
    1100             : static GEN
    1101       37253 : hgmeulerfactorlimit(GEN hgm, GEN t, long p, long d, long flag, long *pe)
    1102             : {
    1103       37253 :   long c = hgmclass(hgm, p, t);
    1104       37256 :   if (c == C_TAME0)
    1105             :   {
    1106          91 :     long v = Q_lvalrem(t, p, &t);
    1107          91 :     return Jordantameexpo(hgm, v, t, p, pe);
    1108             :   }
    1109       37165 :   if (c != C_BAD) return frobpoltrunc(hgm, t, c, p, d, pe);
    1110         301 :   if (flag) { *pe = -1; return gen_0; } else { *pe = 0; return pol_1(0); }
    1111             : }
    1112             : 
    1113             : GEN
    1114         441 : hgmeulerfactor(GEN hgm, GEN t, long p, GEN* pE)
    1115             : {
    1116         441 :   pari_sp av = avma;
    1117             :   long e, B;
    1118             :   GEN P;
    1119         441 :   if (!checkhgm(hgm)) pari_err_TYPE("hgmeulerfactor", hgm);
    1120         441 :   if (!is_rational_t(typ(t))) pari_err_TYPE("hgmeulerfactor",t);
    1121         434 :   if (hgm_get_SWAP(hgm)) t = ginv(t);
    1122         434 :   B = (long)(hgm_get_DEG(hgm) * log(p)) + 1;
    1123         434 :   P = gerepilecopy(av, hgmeulerfactorlimit(hgm, t, p, B, 1, &e));
    1124         420 :   if (pE) *pE = stoi(e);
    1125         420 :   return P;
    1126             : }
    1127             : 
    1128             : /***********************************************************************/
    1129             : /*                        Tame Euler factors                           */
    1130             : /***********************************************************************/
    1131             : /* FIXME: implement properly like RgXn_sqrt */
    1132             : static GEN
    1133          42 : RgXn_sqrtnu(GEN h, long f, long e)
    1134             : {
    1135          42 :   if (f == 1) return h;
    1136           7 :   if (f == 2) return RgXn_sqrt(h, e);
    1137           0 :   return ser2rfrac_i(gsqrtn(RgX_to_ser(h, e + 2), utoipos(f), NULL, 0));
    1138             : }
    1139             : static GEN
    1140          77 : Jordantame(GEN hgm, GEN t0, long m, long p)
    1141             : {
    1142             :   GEN P, T, C, ZP, V;
    1143             :   long d, phim, f, j, c, q, qm, dfp;
    1144             : 
    1145          77 :   if (m == 1)
    1146             :   {
    1147          35 :     long e = hgm_get_WT(hgm) - get_b1(hgm_get_CYCLOE(hgm));
    1148          35 :     return deg1pol_shallow(negi(powuu(p, (e+1) >> 1)), gen_1, 0);
    1149             :   }
    1150          42 :   phim = eulerphiu(m); f = Fl_order(p % m, phim, m);
    1151          42 :   d = phim + 1; /* DEG >= phim >= f */
    1152          42 :   q = upowuu(p, f); qm = (q - 1) / m;
    1153          42 :   V = cgetg(phim + 1, t_VECSMALL);
    1154         133 :   for (j = c = 1; j < m; j++)
    1155          91 :     if (cgcd(j, m) == 1) V[c++] = j * qm;
    1156          42 :   dfp = get_dfp(hgm, p, f);
    1157          42 :   C = hgmCall(hgm, p, f, dfp, V);
    1158          42 :   ZP = zeropadic_shallow(utoipos(p), dfp);
    1159          42 :   T = teich(gadd(t0, ZP)); P = pol_1(0);
    1160         133 :   for (j = 1; j < lg(V); j++)
    1161             :   {
    1162          91 :     GEN Q, c = gmul(gpowgs(T, V[j]), gel(C,j));
    1163          91 :     Q = RgXn_red_shallow(RgX_shift_shallow(RgX_Rg_mul(P, c), f), d);
    1164          91 :     P = RgX_sub(P, Q); /* P * (1 - c x^f) mod X^d */
    1165             :   }
    1166          42 :   return centerlift(RgXn_sqrtnu(P, f, d));
    1167             : }
    1168             : 
    1169             : static GEN
    1170         105 : eulfactameinit(GEN hgm, long v)
    1171             : {
    1172         105 :   GEN C = hgm_get_CYCLOE(hgm);
    1173         105 :   if (!v) pari_err_BUG("hgmeulerfactor [incorrect t in eulfactame]");
    1174         105 :   if (hgm_get_SWAP(hgm)) v = -v;
    1175         105 :   return gel(C, (v < 0)? 2: 1);
    1176             : }
    1177             : static long
    1178          14 : tameexpo(GEN hgm, long v)
    1179             : {
    1180          14 :   GEN W = eulfactameinit(hgm,v);
    1181          14 :   long e, m, l = lg(W);
    1182          42 :   for (m = 1, e = 0; m < l; m++)
    1183          28 :     if (W[m] && v % m == 0) e += eulerphiu(m);
    1184          14 :   return hgm_get_DEG(hgm) - e;
    1185             : }
    1186             : static GEN
    1187          91 : Jordantameexpo(GEN hgm, long v, GEN t0, long p, long *pe)
    1188             : {
    1189          91 :   GEN P = pol_1(0), W = eulfactameinit(hgm, v);
    1190          91 :   long e, m, l = lg(W);
    1191         308 :   for (m = 1, e = 0; m < l; m++)
    1192         217 :     if (W[m] && v % m == 0)
    1193             :     {
    1194          77 :       P = gmul(P, Jordantame(hgm, t0, m, p));
    1195          77 :       e += eulerphiu(m);
    1196             :     }
    1197          91 :   *pe = hgm_get_DEG(hgm) - e; return P;
    1198             : }
    1199             : 
    1200             : /***************************************************************/
    1201             : /*         PART IV.5: Fake wild primes for t                   */
    1202             : /***************************************************************/
    1203             : /* Compute g_q(r)=pi^sq(r)*gauss(r) */
    1204             : static GEN
    1205        3612 : gkgauss(long p, long f, GEN vp, long r, GEN ZP, GEN *sp)
    1206             : {
    1207        3612 :   GEN S = gen_0, P = gen_m1;
    1208        3612 :   long i, qm1 = vp[f+1] - 1;
    1209       13342 :   for (i = 1; i <= f; i++)
    1210             :   {
    1211        9730 :     GEN t = gfrac(sstoQ(r * vp[i], qm1));
    1212        9730 :     S = gadd(S, t);
    1213        9730 :     P = gmul(P, Qp_gamma(gadd(t, ZP)));
    1214             :   }
    1215        3612 :   *sp = gmulsg(p - 1, S); /* t_INT */
    1216        3612 :   return P;
    1217             : }
    1218             : 
    1219             : static GEN
    1220        1428 : hgmG(GEN hgm, long p, long f, GEN vp, long r, GEN ZP)
    1221             : {
    1222        1428 :   GEN S = gen_0, P = gen_1, VPOLGA = hgm_get_VPOLGA(hgm);
    1223             :   long n, c;
    1224        5040 :   for (n = 1; n < lg(VPOLGA); n++)
    1225        3612 :     if ((c = VPOLGA[n]))
    1226             :     {
    1227        3612 :       GEN sq, g = gkgauss(p, f, vp, r*n, ZP, &sq);
    1228        3612 :       S = addii(S, mulsi(c, sq));
    1229        3612 :       P = gmul(P, gpowgs(g, c));
    1230             :     }
    1231             :   /* p - 1 | S */
    1232        1428 :   return gmul(P, powis(utoineg(p), itos(S) / (p - 1)));
    1233             : }
    1234             : 
    1235             : /* multiplicity of -r / (q-1) in beta */
    1236             : static long
    1237        2856 : hgmmulti(GEN B, long q, long r)
    1238             : {
    1239        2856 :   long d = (q-1) / ugcd(r, q-1);
    1240        2856 :   return d >= lg(B)? 0: B[d];
    1241             : }
    1242             : 
    1243             : /* We must have w(M)^r*hgmQ(q,r)=hgmC(q,r)/hgmC(q,0) */
    1244             : static GEN
    1245        1428 : hgmQ(GEN hgm, long p, long f, GEN vp, long r, GEN ZP)
    1246             : {
    1247        1428 :   pari_sp av = avma;
    1248        1428 :   GEN B = gel(hgm_get_CYCLOE(hgm), 2);
    1249        1428 :   long q = vp[f+1], m0 = hgmmulti(B, q, 0), m1 = hgmmulti(B, q, r);
    1250        1428 :   GEN c = powis(utoipos(q), hgm_get_TT(hgm) + m0 - m1);
    1251        1428 :   if (odd(m0)) c = negi(c);
    1252        1428 :   return gerepileupto(av, padic_to_Q(gmul(c, hgmG(hgm, p, f, vp, r, ZP))));
    1253             : }
    1254             : 
    1255             : static GEN
    1256         231 : hgmU(GEN hgm, long p, long f, GEN t, long dfp)
    1257             : {
    1258         231 :   pari_sp av = avma;
    1259         231 :   GEN ZP = zeropadic_shallow(utoipos(p), dfp), vp = upowers_u(p, f, 1);
    1260         231 :   long q = vp[f+1], i;
    1261         231 :   GEN Q = cgetg(q+1, t_POL);
    1262         231 :   Q[1] = evalsigne(1)|evalvarn(0);
    1263        1659 :   for (i = 2; i <= q; i++) gel(Q, i) = hgmQ(hgm, p, f, vp, i, ZP);
    1264         231 :   t = p == 2? gen_1: gmul(hgm_get_MVALUE(hgm), t);
    1265         231 :   return gerepileupto(av, hgmH(Q, p, f, dfp, t));
    1266             : }
    1267             : 
    1268             : /***************************************************************/
    1269             : /*         PART V: Utility programs and main init              */
    1270             : /***************************************************************/
    1271             : static GEN
    1272           7 : cycloE2cyclo(GEN A, GEN B)
    1273           7 : { return mkvec2(count2list(A), count2list(B)); }
    1274             : #if 0
    1275             : GEN
    1276             : hgmalphatocyclo(GEN val, GEN vbe)
    1277             : { GEN C = get_CYCLOE(val,vbe); return cycloE2cyclo(gel(C,1), gel(C,2)); }
    1278             : #endif
    1279             : 
    1280             : static GEN
    1281      843850 : allprims(long n, GEN cache)
    1282             : {
    1283             :   long l, i, c;
    1284             :   GEN z, v;
    1285      843850 :   if (gel(cache,n)) return gel(cache,n);
    1286      418635 :   z = coprimes_zv(n); l = lg(z); v = cgetg(l, t_VEC);
    1287     3036782 :   for (i = c = 1; i < l; i++)
    1288     2618147 :     if (z[i]) gel(v, c++) = mkfracss(i, n);
    1289      418635 :   setlg(v, c); return gel(cache,n) = v;
    1290             : }
    1291             : static GEN
    1292      207956 : zv_to_prims(GEN D, GEN cache)
    1293             : {
    1294      207956 :   long i, l = lg(D);
    1295      207956 :   GEN v = cgetg(l, t_VEC);
    1296     1051806 :   for (i = 1; i < l; i++)
    1297             :   {
    1298      843857 :     if (D[i] <= 0) pari_err_TYPE("hgmcyclotoalpha", D);
    1299      843850 :     gel(v, i) = allprims(D[i], cache);
    1300             :   }
    1301      207949 :   return shallowconcat1(v);
    1302             : }
    1303             : static void
    1304      103978 : hgmcyclotoalpha(GEN *pA, GEN *pB)
    1305             : {
    1306      103978 :   GEN v, D = *pA, E = *pB;
    1307      103978 :   if (typ(D) != t_VECSMALL) D = gtovecsmall(D);
    1308      103978 :   if (typ(E) != t_VECSMALL) E = gtovecsmall(E);
    1309      103978 :   v = const_vec(maxss(vecsmall_max(D), vecsmall_max(E)), NULL);
    1310      103978 :   gel(v,1) = mkvec(gen_0);
    1311      103978 :   *pA = zv_to_prims(D, v);
    1312      103978 :   *pB = zv_to_prims(E, v);
    1313      103971 :   if (lg(*pA) != lg(*pB))
    1314           0 :     pari_err_TYPE("hgminit [incorrect lengths]", mkvec2(D,E));
    1315      103971 : }
    1316             : 
    1317             : static GEN
    1318         231 : hgmalphatogamma(GEN val, GEN vbe) { return get_VPOLGA(get_CYCLOE(val, vbe)); }
    1319             : 
    1320             : static void
    1321          21 : hgmgammatocyclo(GEN VPOLGA, GEN *pD, GEN *pE)
    1322             : {
    1323          21 :   long i, cn, cd, l = lg(VPOLGA);
    1324          21 :   GEN d, n, v = zero_zv(l - 1);
    1325          21 :   if (typ(VPOLGA) != t_VECSMALL) VPOLGA = gtovecsmall(VPOLGA);
    1326         112 :   for (i = 1; i < l; i++)
    1327             :   {
    1328          91 :     long e = VPOLGA[i];
    1329          91 :     if (e)
    1330             :     {
    1331          56 :       GEN D = divisorsu(i);
    1332          56 :       long j, lD = lg(D);
    1333         161 :       for (j = 1; j < lD; j++) v[ D[j] ] += e;
    1334             :     }
    1335             :   }
    1336         112 :   for (i = 1, cn = cd = 0; i < l; i++)
    1337             :   {
    1338          91 :     long e = v[i];
    1339          91 :     if (e > 0) cn += e; else cd -= e;
    1340             :   }
    1341          21 :   if (!cn || !cd) pari_err_TYPE("hgmgammatocyclo", VPOLGA);
    1342          14 :   *pE = n = cgetg(cn+1, t_VECSMALL);
    1343          14 :   *pD = d = cgetg(cd+1, t_VECSMALL);
    1344          98 :   for (i = cn = cd = 1; i < l; i++)
    1345             :   {
    1346          84 :     long e = v[i], j;
    1347         189 :     if (e < 0) for (j = 1; j <= -e; j++) d[cd++] = i;
    1348         126 :     else if (e > 0) for (j = 1; j <= e; j++) n[cn++] = i;
    1349             :   }
    1350          14 : }
    1351             : 
    1352             : static void
    1353          21 : hgmgammatoalpha(GEN VPOLGA, GEN *pA, GEN *pB)
    1354          21 : { hgmgammatocyclo(VPOLGA, pA, pB); hgmcyclotoalpha(pA, pB); }
    1355             : 
    1356             : /* A and B sorted */
    1357             : static long
    1358      554932 : zv_intersect(GEN A, GEN B)
    1359             : {
    1360      554932 :   long a = 1, b = 1, lA = lg(A), lB = lg(B);
    1361     1801464 :   while (a < lA && b < lB)
    1362             :   {
    1363     1697794 :     if      (A[a] < B[b]) a++;
    1364     1379231 :     else if (A[a] > B[b]) b++; else return 1;
    1365             :   }
    1366      103670 :   return 0;
    1367             : }
    1368             : static void
    1369         231 : remove_intersect(GEN *pA, GEN *pB)
    1370             : {
    1371         231 :   GEN V, W, A = *pA, B = *pB;
    1372         231 :   long a = 1, b = 1, v = 1, w = 1, lA, lB;
    1373         231 :   *pA = V = cgetg_copy(A, &lA);
    1374         231 :   *pB = W = cgetg_copy(B, &lB);
    1375        1267 :   while (a < lA && b < lB)
    1376             :   {
    1377        1036 :     int s = gcmp(gel(A,a), gel(B,b));
    1378        1036 :     if   (s < 0) gel(V,v++) = gel(A,a++);
    1379         707 :     else if (s > 0) gel(W,w++) = gel(B,b++);
    1380         273 :     else { a++; b++; }
    1381             :   }
    1382         364 :   while (a < lA) gel(V,v++) = gel(A,a++);
    1383         259 :   while (b < lB) gel(W,w++) = gel(B,b++);
    1384         231 :   setlg(V,v); setlg(W,w);
    1385         231 : }
    1386             : 
    1387             : /* al, be normalized, sorted, unknown intersection */
    1388             : static GEN
    1389         231 : albe2u(GEN al, GEN be)
    1390             : {
    1391             :   GEN ga;
    1392         231 :   remove_intersect(&al, &be);
    1393         231 :   ga = hgmalphatogamma(al, be);
    1394         231 :   return F2v_factorback(ga);
    1395             : }
    1396             : static GEN
    1397         462 : get_u(GEN al, GEN be, GEN CYCLOE, GEN VPOLGA, long DEG, long WT)
    1398             : {
    1399         462 :   GEN u, u0 = F2v_factorback(VPOLGA);
    1400         462 :   long b1 = get_b1(CYCLOE);
    1401         462 :   if (odd(DEG))
    1402             :   {
    1403         147 :     be = QV_normalize(RgV_addhalf(be));
    1404         147 :     u = albe2u(al, be);
    1405         147 :     if ((DEG & 3) == 3) u = negi(u);
    1406             :   }
    1407         315 :   else if (odd(WT)) { u = u0; if ((b1 & 3) == 2) u = negi(u); }
    1408             :   else
    1409             :   {
    1410          84 :     al = QV_normalize(RgV_addhalf(al));
    1411          84 :     u = shifti(albe2u(al, be), 1);
    1412          84 :     if (((DEG + b1) & 3) == 1) u = negi(u);
    1413             :   }
    1414         462 :   u0 = shifti(u0, 1); if ((b1 & 3) == 3) u0 = negi(u0);
    1415         462 :   return mkvec2(coredisc(u), u0);
    1416             : }
    1417             : 
    1418             : static long
    1419          35 : zv_sumeuler(GEN v)
    1420             : {
    1421          35 :   long i, l = lg(v);
    1422          35 :   GEN s = gen_0;
    1423         133 :   for (i = 1; i < l; i++)
    1424             :   {
    1425          98 :     if (v[i] <= 0) pari_err_TYPE("hgminit", v);
    1426          98 :     s = addiu(s, eulerphiu(v[i]));
    1427             :   }
    1428          35 :   return itou(s);
    1429             : }
    1430             : 
    1431             : /* (a, b): format (1) if denominator, format (2) if no denominator,
    1432             :  * format (3) if b not vector. */
    1433             : static GEN
    1434         511 : hgminit_i(GEN a, GEN b)
    1435             : {
    1436         511 :   long ta = typ(a), l = lg(a);
    1437         511 :   if (ta != t_VEC && ta != t_VECSMALL) pari_err_TYPE("hgminit", a);
    1438         504 :   if (!b)
    1439             :   {
    1440          84 :     if (l == 1) initab(a, a); /* error */
    1441          77 :     if (ta == t_VECSMALL || RgV_is_ZV(a))
    1442          49 :     {
    1443             :       long i;
    1444          56 :       if (ta != t_VECSMALL) a = vec_to_vecsmall(a);
    1445         161 :       for (i = 1; i < l; i++)
    1446         126 :         if (a[i] <= 0) break;
    1447          56 :       if (i != l)
    1448          21 :         hgmgammatoalpha(a, &a, &b); /* gamma */
    1449             :       else
    1450             :       { /* cyclo */
    1451          35 :         b = const_vecsmall(zv_sumeuler(a), 1);
    1452          35 :         hgmcyclotoalpha(&a, &b);
    1453             :       }
    1454             :     }
    1455             :     else /* alpha */
    1456          21 :       b = zerovec(l - 1);
    1457             :   }
    1458             :   else
    1459             :   {
    1460         420 :     if (typ(b) != ta) pari_err_TYPE("hgminit", b);
    1461         413 :     if (l > 1 && (ta == t_VECSMALL || (RgV_is_ZV(a) && RgV_is_ZV(b))))
    1462         259 :       hgmcyclotoalpha(&a, &b); /* cyclo */
    1463             :   }
    1464         476 :   return initab(a, b);
    1465             : }
    1466             : GEN
    1467         511 : hgminit(GEN val, GEN vbe)
    1468         511 : { pari_sp av = avma; return gerepilecopy(av, hgminit_i(val, vbe)); }
    1469             : 
    1470             : GEN
    1471          14 : hgmalpha(GEN hgm)
    1472             : {
    1473          14 :   GEN al = hgm_get_VAL(hgm), be = hgm_get_VBE(hgm);
    1474          14 :   if (hgm_get_SWAP(hgm)) swap(al, be);
    1475          14 :   retmkvec2(gcopy(al), gcopy(be));
    1476             : }
    1477             : GEN
    1478          14 : hgmgamma(GEN hgm)
    1479             : {
    1480          14 :   pari_sp av = avma;
    1481             :   GEN g;
    1482          14 :   if (!checkhgm(hgm)) pari_err_TYPE("hgmgamma", hgm);
    1483          14 :   g = hgm_get_VPOLGA(hgm);
    1484          14 :   if (hgm_get_SWAP(hgm)) g = zv_neg(g);
    1485          14 :   return gerepilecopy(av, g);
    1486             : }
    1487             : GEN
    1488           7 : hgmcyclo(GEN hgm)
    1489             : {
    1490           7 :   pari_sp av = avma;
    1491             :   GEN A, B, C;
    1492           7 :   if (!checkhgm(hgm)) pari_err_TYPE("hgmcyclo", hgm);
    1493           7 :   C = hgm_get_CYCLOE(hgm); A = gel(C,1); B = gel(C,2);
    1494           7 :   if (hgm_get_SWAP(hgm)) swap(A, B);
    1495           7 :   return gerepilecopy(av, cycloE2cyclo(A, B));
    1496             : }
    1497             : 
    1498             : GEN
    1499          21 : hgmtwist(GEN hgm)
    1500             : {
    1501          21 :   pari_sp av = avma;
    1502             :   GEN val, vbe;
    1503          21 :   if (!checkhgm(hgm)) pari_err_TYPE("hgmtwist", hgm);
    1504          21 :   val = hgm_get_VAL(hgm);
    1505          21 :   vbe = hgm_get_VBE(hgm);
    1506          21 :   if (hgm_get_SWAP(hgm)) swap(val, vbe);
    1507          21 :   val = sort(RgV_addhalf(val));
    1508          21 :   vbe = sort(RgV_addhalf(vbe));
    1509          21 :   return gerepilecopy(av, initab(val, vbe));
    1510             : }
    1511             : 
    1512             : GEN
    1513         161 : hgmparams(GEN hgm)
    1514             : {
    1515         161 :   pari_sp av = avma;
    1516             :   GEN H, M;
    1517             :   long TT, DEG, WT;
    1518         161 :   if (!checkhgm(hgm)) pari_err_TYPE("hgmparams", hgm);
    1519         161 :   H = zx_to_ZX(hgm_get_HODGE(hgm));
    1520         161 :   TT = hgm_get_TT(hgm); DEG = hgm_get_DEG(hgm);
    1521         161 :   WT = hgm_get_WT(hgm); M = hgm_get_MVALUE(hgm);
    1522         161 :   return gerepilecopy(av, mkvec4(utoipos(DEG), utoi(WT),
    1523             :                                  mkvec2(H,stoi(TT)), M));
    1524             : }
    1525             : 
    1526             : /* symmetry at 1? */
    1527             : long
    1528           7 : hgmissymmetrical(GEN hgm)
    1529             : {
    1530             :   GEN A, B, C;
    1531             :   long a, i, j, lA, lB;
    1532           7 :   if (!checkhgm(hgm)) pari_err_TYPE("hgmissymmetrical", hgm);
    1533           7 :   C = hgm_get_CYCLOE(hgm);
    1534           7 :   if (odd(hgm_get_DEG(hgm))) return 0;
    1535           7 :   A = gel(C,1); lA = lg(A);
    1536           7 :   B = gel(C,2); lB = lg(B);
    1537          21 :   for (i = 1; i < lA; i++) if ((a = A[i]))
    1538             :   {
    1539           7 :     switch(i & 3)
    1540             :     { /* polcyclo(i)[-X] = polcyclo(j) */
    1541           0 :       case 0: j = i; break;
    1542           7 :       case 2: j = i >> 1; break;
    1543           0 :       default: j = i << 1; break;
    1544             :     }
    1545           7 :     if (j >= lB || B[j] != a) return 0;
    1546             :   }
    1547           7 :   return 1;
    1548             : }
    1549             : 
    1550             : /***************************************************************/
    1551             : /*         PART VI: Experimental euler factors                 */
    1552             : /***************************************************************/
    1553             : /* FIXME: one prime at a time */
    1554             : static GEN
    1555      227731 : hgmmodif(GEN an, GEN PPol)
    1556             : {
    1557      227731 :   pari_sp av = avma;
    1558      227731 :   long i, L = lg(an), lP = lg(PPol);
    1559             : 
    1560      674401 :   for (i = 1; i < lP; i++)
    1561             :   {
    1562      446670 :     GEN E = gel(PPol, i), pol = gel(E, 2);
    1563      446670 :     long d = degpol(pol);
    1564      446670 :     if (d)
    1565             :     {
    1566      433629 :       GEN v = vec_ei(L, 1);
    1567      433629 :       long j, p = itos(gel(E, 1)), pj = p;
    1568     1384271 :       for (j = 1; j <= d && pj <= L; j++, pj *= p)
    1569      950642 :         gel(v, pj) = RgX_coeff(pol, j);
    1570      433629 :       an = dirdiv(an, v);
    1571             :     }
    1572             :   }
    1573      227731 :   return gerepileupto(av, an);
    1574             : }
    1575             : 
    1576             : /***************************************************************/
    1577             : /*             PART VII: Make tables of HGMs                   */
    1578             : /***************************************************************/
    1579             : 
    1580             : static int
    1581         378 : ok_part(GEN v, GEN W)
    1582             : {
    1583         378 :   long l = lg(v), j;
    1584        1617 :   for (j = 1; j < l; j++)
    1585        1435 :     if (!gel(W,v[j])) return 0;
    1586         182 :   return 1;
    1587             : }
    1588             : static GEN
    1589          21 : mkphi()
    1590             : {
    1591          21 :   GEN W = const_vec(20, NULL);
    1592          21 :   gel(W,1) = mkvecsmall2(1,2);
    1593          21 :   gel(W,2) = mkvecsmall3(3,4,6);
    1594          21 :   gel(W,4) = mkvecsmall4(5,8,10,12);
    1595          21 :   gel(W,6) = mkvecsmall4(7,9,14,18);
    1596          21 :   gel(W,8) = mkvecsmall5(15,16,20,24,30);
    1597          21 :   gel(W,10)= mkvecsmall2(11,22);
    1598          21 :   gel(W,12)= mkvecsmalln(6, 13L,21L,26L,28L,36L,42L);
    1599          21 :   gel(W,16)= mkvecsmalln(6, 17L,32L,34L,40L,48L,60L);
    1600          21 :   gel(W,18)= mkvecsmall4(19,27,38,54);
    1601          21 :   gel(W,20)= mkvecsmall5(25,33,44,50,66);
    1602          21 :   return W;
    1603             : }
    1604             : 
    1605             : static GEN
    1606         182 : mkal(GEN p, GEN W)
    1607             : {
    1608             :   GEN res, V;
    1609         182 :   long l = lg(p), lV, i, j;
    1610         182 :   V = cgetg(l, t_VECSMALL);
    1611         987 :   for (i = 1; i < l; i++) V[i] = lg(gel(W, p[i])) - 1;
    1612         182 :   V = cyc2elts(V); lV = lg(V);
    1613         182 :   res = cgetg(1, t_VEC);
    1614       20769 :   for (j = 1; j < lV; j++)
    1615             :   {
    1616       20587 :     GEN t = cgetg(l, t_VECSMALL), v = gel(V,j);
    1617      157185 :     for (i = 1; i < l; i++) t[i] = umael2(W, p[i], v[i]+1);
    1618       20587 :     vecsmall_sort(t); if (!RgV_isin(res, t)) res = vec_append(res, t);
    1619             :   }
    1620         182 :   return res;
    1621             : }
    1622             : static GEN
    1623          21 : mkallal(long n)
    1624             : {
    1625          21 :   GEN W = mkphi(), p = partitions(n, NULL, NULL);
    1626          21 :   long i, c, l = lg(p);
    1627         399 :   for (i = c = 1; i < l; i++)
    1628         378 :     if (ok_part(gel(p,i), W)) gel(p, c++) = mkal(gel(p,i), W);
    1629          21 :   setlg(p, c); return shallowconcat1(p);
    1630             : }
    1631             : 
    1632             : static GEN
    1633          21 : mkalbe(long n)
    1634             : {
    1635          21 :   GEN w, L = mkallal(n);
    1636          21 :   long i, j, c, l = lg(L);
    1637          21 :   w = cgetg(l * (l / 2) + 1, t_VEC);
    1638        3941 :   for (i = c = 1; i < l; i++)
    1639             :   {
    1640        3920 :     GEN A = gel(L, i);
    1641        3920 :     long a = A[lg(A)-1];
    1642      558852 :     for (j = i + 1; j < l; j++)
    1643             :     {
    1644      554932 :       GEN B = gel(L, j);
    1645      554932 :       if (!zv_intersect(A, B))
    1646      179081 :         gel(w,c++) = (A[1]==1 || (B[1]!=1 && a > B[lg(B)-1]))?
    1647      179081 :                      mkvec2(B,A): mkvec2(A,B);
    1648             :     }
    1649             :   }
    1650          21 :   setlg(w,c); return w;
    1651             : }
    1652             : 
    1653             : static long
    1654      103670 : cyclowt(GEN a, GEN b)
    1655             : {
    1656      103670 :   pari_sp av = avma;
    1657             :   long TT;
    1658      103670 :   hgmcyclotoalpha(&a, &b);
    1659      103670 :   return gc_long(av, degpol(hodge(a, b, &TT)));
    1660             : }
    1661             : 
    1662             : GEN
    1663          21 : hgmbydegree(long n)
    1664             : {
    1665          21 :   pari_sp av = avma;
    1666          21 :   GEN w = cgetg(n+1, t_VEC), c = const_vecsmall(n, 1), v = mkalbe(n);
    1667          21 :   long i, l = lg(v);
    1668         154 :   for (i = 1; i <= n; i++) gel(w,i) = cgetg(l,t_VEC);
    1669      103691 :   for (i = 1; i < l; i++)
    1670             :   {
    1671      103670 :     GEN z = gel(v,i);
    1672      103670 :     long k = cyclowt(gel(z,1), gel(z,2)) + 1;
    1673      103670 :     gmael(w, k, c[k]++) = z;
    1674             :   }
    1675         154 :   for (i = 1; i <= n; i++) setlg(gel(w,i), c[i]);
    1676          21 :   return gerepilecopy(av, w);
    1677             : }
    1678             : 
    1679             : /***************************************************************/
    1680             : /*             PART VIII: L-function data                      */
    1681             : /***************************************************************/
    1682             : /* BAD sorted t_VECSMALL */
    1683             : static GEN
    1684         217 : removebad(GEN v, GEN BAD)
    1685             : {
    1686         217 :   long i, c, l = lg(v);
    1687         217 :   GEN w = cgetg(l, t_VECSMALL);
    1688         434 :   for (i = c = 1; i < lg(v); i++)
    1689         217 :     if (!zv_search(BAD, v[i])) w[c++] = v[i];
    1690         217 :   setlg(w, c); return w;
    1691             : }
    1692             : static GEN
    1693         609 : primedivisors(GEN t)
    1694         609 : { GEN fa = absZ_factor(t); return ZV_to_zv(gel(fa,1)); }
    1695             : 
    1696             : static GEN
    1697         385 : gacfac(long p, long m, long c)
    1698             : {
    1699         385 :   long i, l = p + m + 1;
    1700         385 :   GEN v = cgetg(l, t_VECSMALL);
    1701         728 :   for (i = 1; i <= p; i++) v[i] = -c;
    1702         658 :   for (     ; i <  l; i++) v[i] = 1 - c;
    1703         385 :   return v;
    1704             : }
    1705             : 
    1706             : static GEN
    1707         203 : hgmfindvga(GEN hgm, GEN t)
    1708             : {
    1709         203 :   GEN v, HODGE = hgm_get_HODGE(hgm);
    1710         203 :   long WT = degpol(HODGE), WT2 = (WT - 1) >> 1, i, c, h, fl = gequal1(t);
    1711         203 :   v = cgetg(WT + 2, t_VEC);
    1712         413 :   for (i = 0, c = 1; i <= WT2; i++)
    1713             :   {
    1714         210 :     if ((h = HODGE[i+2]))
    1715             :     {
    1716         210 :       if (fl && 2*i == WT - 1) h--;
    1717         210 :       if (h) gel(v, c++) = gacfac(h, h, i);
    1718             :     }
    1719             :   }
    1720         203 :   if (!odd(WT))
    1721             :   {
    1722         182 :     long h = HODGE[WT2+3], hp = h >> 1, hm;
    1723         182 :     if (!fl)
    1724             :     {
    1725         182 :       hm = hp;
    1726         182 :       if (odd(h))
    1727         154 :       { if (gsigne(t) >= 0 && gcmpgs(t, 1) <= 0) hm++; else hp++; }
    1728             :       else
    1729          28 :       { if (gcmpgs(t, 1) > 0) { hp++; hm--; } }
    1730         182 :       if (odd(hgm_get_TT(hgm) + WT2 + 1)) lswap(hp, hm);
    1731             :     }
    1732             :     else
    1733             :     {
    1734           0 :       h--; hm = h - hp;
    1735           0 :       if (odd(h) && odd(hgm_get_TT(hgm) + WT2 + 1)) lswap(hp, hm);
    1736             :     }
    1737         182 :     if (hp || hm) gel(v, c++) = gacfac(hp, hm, WT2 + 1);
    1738             :   }
    1739         203 :   if (c == 1) return cgetg(1, t_VECSMALL);
    1740         203 :   setlg(v, c); v = shallowconcat1(v);
    1741         203 :   vecsmall_sort(v); return v;
    1742             : }
    1743             : 
    1744             : /* Return [VGA, k, BAD, COND] where VGA as in lfun, k non motivic
    1745             :  * weight (s -> k - s), BAD is the list of wild primes and Euler factors,
    1746             :  * COND is the tame part of the conductor */
    1747             : static GEN
    1748         203 : hgmlfuninfty(GEN hgm, GEN t)
    1749             : {
    1750         203 :   pari_sp av = avma;
    1751         203 :   GEN VGA = hgmfindvga(hgm, t), T0, T1, v;
    1752         203 :   GEN COND, t1 = gsubgs(t, 1), BAD = primedivisors(hgm_get_BAD(hgm));
    1753         203 :   long k = hgm_get_WT(hgm) + 1, i, j, l;
    1754             : 
    1755         203 :   if (isintzero(t1)) T1 = cgetg(1, t_VECSMALL);
    1756             :   else
    1757             :   {
    1758         196 :     GEN fa = absZ_factor(numer_i(t1)), P = gel(fa,1), E = gel(fa,2);
    1759         196 :     if (!odd(k)) T1 = removebad(ZV_to_zv(P), BAD);
    1760             :     else
    1761             :     {
    1762         182 :       long j, lP = lg(P);
    1763         182 :       T1 = cgetg(lP, t_VECSMALL);
    1764         350 :       for (i = j = 1; i < lP; i++)
    1765             :       {
    1766         168 :         long p = itos(gel(P,i));
    1767         168 :         if (mpodd(gel(E,i)) && !zv_search(BAD, p)) T1[j++] = p;
    1768             :       }
    1769         182 :       setlg(T1,j);
    1770             :     }
    1771             :   }
    1772         203 :   COND = zv_prod_Z(T1);
    1773         203 :   T0 = removebad(gconcat(primedivisors(numer_i(t)),
    1774             :                          primedivisors(denom_i(t))), BAD);
    1775         217 :   for (i = 1; i < lg(T0); i++)
    1776             :   {
    1777          14 :     long p = T0[i], e = tameexpo(hgm, Q_lval(t, p));
    1778          14 :     COND = mulii(powuu(p, e), COND);
    1779             :   }
    1780         203 :   l = lg(BAD); v = cgetg(l, t_VEC);
    1781         490 :   for (i = j = 1; i < l; i++) /* FIXME: precompute wild Euler factors */
    1782         287 :     if (hgmclass(hgm, BAD[i], t) == C_BAD)
    1783         252 :       gel(v,j++) = mkvec2(utoipos(BAD[i]), pol_1(0));
    1784         203 :   setlg(v, j); return gerepilecopy(av, mkvec4(VGA, utoi(k), v, COND));
    1785             : }
    1786             : 
    1787             : /***************************************************************/
    1788             : /*             PART IX: GUESS OF WILD PRIMES                   */
    1789             : /***************************************************************/
    1790             : static GEN
    1791        1708 : vecmellin(GEN L, GEN K, GEN t, GEN td, GEN vroots, long bitprec)
    1792             : {
    1793        1708 :   long i, N = lfunthetacost(L, t, 0, bitprec);
    1794        1708 :   GEN v = cgetg(N+1, t_VEC);
    1795      427231 :   for (i = 1; i <= N; i++)
    1796      425523 :     gel(v,i) = gammamellininvrt(K, gmul(td, gel(vroots,i)), bitprec);
    1797        1708 :   return v;
    1798             : }
    1799             : 
    1800             : /* List of Weil polynomials for prime p of degree d and weight w */
    1801             : static GEN
    1802        1610 : listweil_i(ulong d, ulong p, ulong w)
    1803             : {
    1804             :   GEN V;
    1805        1610 :   if (d == 0) return mkvec(pol_1(0));
    1806        1225 :   if (odd(d))
    1807             :   {
    1808             :     GEN P;
    1809         644 :     if (odd(w)) return cgetg(1, t_VEC);
    1810         420 :     V = listweil_i(d - 1, p, w);
    1811         420 :     P = monomial(powuu(p, w / 2), 1, 0);
    1812         420 :     return shallowconcat(gmul(gsubsg(1, P), V), gmul(gaddsg(1, P), V));
    1813             :   }
    1814         581 :   if (d == 2)
    1815             :   {
    1816         553 :     long q = upowuu(p, w), s = usqrt(4 * q), i;
    1817         553 :     GEN Q = utoi(q);
    1818         553 :     V = cgetg(2*s + 3, t_VEC);
    1819        4844 :     for (i = 1; i <= 2*s + 1; i++) gel(V,i) = mkpoln(3, Q, stoi(1+s-i), gen_1);
    1820         553 :     gel(V, 2*s + 2) = mkpoln(3, negi(Q), gen_0, gen_1);
    1821         553 :     return V;
    1822             :   }
    1823          28 :   if (d == 4)
    1824             :   {
    1825          28 :     long q = upowuu(p, w), s = usqrt(16 * q), a, c, is2 = usqrt(4 * q);
    1826          28 :     double s2 = 2 * sqrt((double)q);
    1827          28 :     GEN Q2 = sqru(q), W, A, mA, tmp, Q;
    1828          28 :     V = cgetg(s + 3, t_VEC);
    1829         252 :     for (a = 0; a <= s; a++)
    1830             :     {
    1831         224 :       long b, m = ceil(a * s2) - 2 * q, M = ((a * a) >> 2) + 2 * q;
    1832         224 :       GEN AQ = stoi(a * q), mAQ = stoi(-a * q);
    1833         224 :       A = stoi(a); mA = stoi(-a);
    1834         224 :       W = cgetg(2*(M - m) + 3, t_VEC);
    1835        1876 :       for (b = m, c = 1; b <= M; b++)
    1836             :       {
    1837        1652 :         if (a) gel(W, c++) = mkpoln(5, Q2, mAQ, stoi(b), mA, gen_1);
    1838        1652 :         gel(W, c++) = mkpoln(5, Q2, AQ, stoi(b), A, gen_1);
    1839             :       }
    1840         224 :       setlg(W, c); gel(V, a + 1) = W;
    1841             :     }
    1842          28 :     W = cgetg(2 * is2 + 2, t_VEC);
    1843          28 :     tmp = mkpoln(3, stoi(-q), gen_0, gen_1);
    1844          28 :     Q = utoipos(q);
    1845         147 :     for (a = 0, c = 1; a <= is2; a++)
    1846             :     {
    1847         119 :       A = stoi(a); mA = stoi(-a);
    1848         119 :       if (a) gel(W, c++) = gmul(tmp, mkpoln(3, Q, mA, gen_1));
    1849         119 :       gel(W, c++) = gmul(tmp, mkpoln(3, Q, A, gen_1));
    1850             :     }
    1851          28 :     setlg(W, c); gel(V, s + 2) = W;
    1852          28 :     return shallowconcat1(V);
    1853             :   }
    1854           0 :   pari_err_IMPL("d > 5 in listweil");
    1855             :   return NULL; /* LCOV_EXCL_LINE */
    1856             : }
    1857             : 
    1858             : /* Same, weight <= w, by decreasing weight */
    1859             : static GEN
    1860         490 : listweilallw_i(ulong d, ulong p, ulong w)
    1861             : {
    1862         490 :   long i = d? w: 0;
    1863         490 :   GEN V = cgetg(i+2, t_VEC);
    1864        1680 :   for (; i >= 0; i--) gel(V,i+1) = listweil_i(d, p, i);
    1865         490 :   return shallowconcat1(V);
    1866             : }
    1867             : 
    1868             : static long
    1869      127344 : sumdeg(GEN W, GEN M)
    1870             : {
    1871      127344 :   long i, l = lg(M), s = 0;
    1872      375851 :   for (i = 1; i < l; i++) s += degpol(gmael(W,i,M[i]+1));
    1873      127344 :   return s;
    1874             : }
    1875             : 
    1876             : static GEN
    1877          35 : mul(GEN a, GEN b, ulong lim)
    1878             : {
    1879          35 :   long na = lg(a)-1, nb = lg(b)-1, i, j, n;
    1880          35 :   GEN c = cgetg(na * nb + 1, t_VECSMALL);
    1881         378 :   for (i = n = 1; i <= na; i++)
    1882        2597 :     for (j = 1; j <= nb; j++)
    1883             :     {
    1884        2254 :       ulong m = umuluu_or_0(a[i], b[j]);
    1885        2254 :       if (m && m <= lim) c[n++] = m;
    1886             :     }
    1887          35 :   setlg(c, n); return c;
    1888             : }
    1889             : static GEN
    1890          98 : listcond(GEN BAD, GEN achi, ulong min, ulong max)
    1891             : {
    1892          98 :   long i, j, l = lg(BAD);
    1893          98 :   GEN P = cgetg(l, t_VEC), z;
    1894         231 :   for (i = 1; i < l; i++)
    1895             :   {
    1896         133 :     long p = BAD[i], v = achi[i];
    1897         133 :     gel(P,i) = upowers_u(p, ulogint(max, p) - v, upowuu(p, v));
    1898             :   }
    1899         133 :   z = gel(P,1); for (i = 2; i < l; i++) z = mul(z, gel(P,i), max);
    1900          98 :   if (min > 1)
    1901             :   {
    1902           0 :     l = lg(z);
    1903           0 :     for (i = j = 1; i < l; i++)
    1904           0 :       if ((ulong)z[i] >= min) z[j++] = z[i];
    1905           0 :     setlg(z, j);
    1906             :   }
    1907          98 :   vecsmall_sort(z); return z;
    1908             : }
    1909             : 
    1910             : /* Artificial lfundiv by zeta(s - (k - 1) / 2). */
    1911             : static GEN
    1912          35 : lfundivraw(GEN L)
    1913             : {
    1914          35 :   long k = itos(ldata_get_k(L));
    1915          35 :   L = shallowcopy(L);
    1916          35 :   gel(L, 3) = lfunvgasub(ldata_get_gammavec(L), mkvec(stoi(-(k - 1) >> 2)));
    1917          35 :   setlg(L, 7); return L;
    1918             : }
    1919             : 
    1920             : /* Compute forvec vectors from v, sorted by increasing total degree */
    1921             : static GEN
    1922          98 : forvecsort(GEN vF, GEN v)
    1923             : {
    1924          98 :   GEN w, E = cyc2elts(v);
    1925          98 :   long i, l = lg(E);
    1926          98 :   w = cgetg(l, t_VECSMALL);
    1927      127442 :   for (i = 1; i < l; i++) w[i] = sumdeg(vF, gel(E, i));
    1928          98 :   return vecpermute(E, vecsmall_indexsort(w)); /* by increasing total degree */
    1929             : }
    1930             : 
    1931             : static GEN
    1932          98 : lfunhgmwild(GEN L, GEN H, GEN t, GEN BAD, long pole, GEN hint, long bitprec)
    1933             : {
    1934             :   GEN v, K, t0, t0r, t0ir, t0i, t0k, N0, vM, vD, val, PPOL, vF, achi;
    1935          98 :   long d, lM, iN, iM, i, k, k2, prec = nbits2prec(bitprec), lB = lg(BAD);
    1936             :   long BADprod, limdeg;
    1937          98 :   ulong minN = 1, maxN = 2048;
    1938             : 
    1939          98 :   v = cgetg(lB, t_VECSMALL); PPOL = cgetg(lB, t_VEC);
    1940         231 :   for (i = 1; i < lB; i++)
    1941             :   {
    1942         133 :     v[i] = itou( gmael(BAD,i,1) );
    1943         133 :     gel(PPOL,i) = shallowcopy(gel(BAD,i));
    1944             :   }
    1945          98 :   BAD = v;
    1946          98 :   BADprod = zv_prod(BAD);
    1947          98 :   achi = get_achi(H, t, BAD);
    1948          98 :   if (pole) L = lfundivraw(L);
    1949          98 :   limdeg = d = ldata_get_degree(L);
    1950          98 :   N0 = ldata_get_conductor(L);
    1951          98 :   if (hint) switch(typ(hint))
    1952             :   {
    1953             :     long l;
    1954           0 :     case t_INT:
    1955           0 :       limdeg = itos(hint);
    1956           0 :       if (limdeg < 0 || limdeg > d) pari_err_TYPE("lfunhgm [bad hint]", hint);
    1957           0 :       break;
    1958           0 :     case t_VEC:
    1959           0 :       l = lg(hint);
    1960           0 :       if (l > 1 && l < 4)
    1961             :       {
    1962           0 :         GEN t = gel(hint,1), r;
    1963           0 :         if (typ(t) != t_INT || signe(t) <= 0)
    1964           0 :           pari_err_TYPE("lfunhgm [bad hint]", hint);
    1965           0 :         t = dvmdii(t, N0, &r);
    1966           0 :         if (r != gen_0)
    1967           0 :           pari_err_TYPE("lfunhgm [bad hint]", hint);
    1968           0 :         minN = maxN = itou(t);
    1969           0 :         if (l == 3)
    1970             :         {
    1971           0 :           t = gel(hint,2);
    1972           0 :           if (typ(t) != t_INT || signe(t) < 0 || cmpis(t, d) > 0)
    1973           0 :             pari_err_TYPE("lfunhgm [bad hint]", hint);
    1974           0 :           limdeg = itos(t);
    1975             :         }
    1976             :       }
    1977             :   }
    1978          98 :   k = itos(ldata_get_k(L)); k2 = (k-1) >> 1;
    1979          98 :   K = gammamellininvinit(ldata_get_gammavec(L), 0, bitprec + 32);
    1980          98 :   t0 = sstoQ(11, 10); t0i = ginv(t0); t0k = gpowgs(t0, k);
    1981          98 :   t0r = gpow(t0, sstoQ(2,d), prec); t0ir = ginv(t0r);
    1982             :   /* vF[i] list of Weil polynomials F for p = BAD[i],
    1983             :    * F = prod_j Fj, |reciprocal root of Fj|^2 = p^(w + 1 - nj)
    1984             :    * 2v_p(lc(F)) = deg F * (w + 1) - sum_j deg Fj * nj; */
    1985          98 :   vF = cgetg(lB, t_VEC);
    1986          98 :   vD = cgetg(lB, t_VEC); /* vD[k][l] = sum_j deg Fj * nj for F = vF[k][l] */
    1987          98 :   v = cgetg(lB, t_VECSMALL);
    1988         231 :   for (i = 1; i < lB; i++)
    1989             :   {
    1990         133 :     GEN W = cgetg(limdeg+2, t_VEC), D;
    1991         133 :     long p = BAD[i], j, lW;
    1992         623 :     for (j = 0; j <= limdeg; j++) gel(W,j+1) = listweilallw_i(j, p, k-1);
    1993         133 :     gel(vF, i) = W = shallowconcat1(W);
    1994         133 :     lW = lg(W); v[i] = lW-1;
    1995         133 :     gel(vD,i) = D = cgetg(lW, t_VEC);
    1996       10136 :     for (j = 1; j < lW; j++)
    1997             :     {
    1998       10003 :       GEN F = gel(W,j);
    1999       10003 :       D[j] = degpol(F) * k - 2 * Z_lval(leading_coeff(F), p);
    2000             :     }
    2001             :   }
    2002          98 :   vM = forvecsort(vF, v); lM = lg(vM);
    2003          98 :   if (DEBUGLEVEL) { err_printf(" lM = %ld ", lM); err_flush(); }
    2004          98 :   L = shallowcopy(L);
    2005          98 :   val = cgetg(lB, t_VECSMALL);
    2006             :   for(;;)
    2007           0 :   {
    2008             :     GEN z, vroots, an0, vN;
    2009             :     long lN, lim;
    2010             : 
    2011          98 :     z = listcond(BAD, achi, minN, maxN);
    2012          98 :     if (maxN == minN) /* from hint */
    2013             :     {
    2014           0 :       minN = 1; /* in case it fails */
    2015           0 :       maxN--;
    2016             :     }
    2017             :     else
    2018             :     {
    2019          98 :       minN = maxN+1;
    2020          98 :       maxN <<= 2;
    2021             :     }
    2022          98 :     lN = lg(z);
    2023          98 :     if (lN == 1) continue;
    2024          98 :     vN = cgetg(lN, t_VEC);
    2025        1785 :     for (i = 1; i < lN; i++) gel(vN, i) = muliu(N0, z[i]);
    2026          98 :     gel(L, 5) = gel(vN, lN - 1);
    2027          98 :     if (DEBUGLEVEL) err_printf("\nmaxN = %ld\n", itos(gel(L,5)));
    2028          98 :     lim = lfunthetacost(L, t0i, 0, bitprec);
    2029          98 :     vroots = dirpowers(lim, sstoQ(2, d), prec + EXTRAPRECWORD);
    2030          98 :     an0 = hgmcoefs(H, t, lim);
    2031          98 :     if (pole)
    2032             :     {
    2033          35 :       GEN w = vecpowuu(lim, k2);
    2034       44450 :       for (i = 1; i <= lim; i++)
    2035       44415 :         if (cgcd(i, BADprod) > 1) gel(w, i) = gen_0;
    2036          35 :       an0 = dirdiv(an0, w);
    2037             :     }
    2038         854 :     for (iN = 1; iN < lN; iN++)
    2039             :     {
    2040         854 :       pari_sp av0 = avma, av2;
    2041         854 :       GEN vecK, vecKi, N = gel(vN,iN), isqN = gpow(N, sstoQ(-1,d), prec);
    2042         854 :       if (DEBUGLEVEL) { err_printf(" %ld ", itos(N)); err_flush(); }
    2043         854 :       gel(L,5) = N;
    2044         854 :       vecK = vecmellin(L, K, t0, gmul(t0r,  isqN), vroots, bitprec);
    2045         854 :       vecKi= vecmellin(L, K, t0i,gmul(t0ir, isqN), vroots, bitprec);
    2046        2310 :       for (i = 1; i < lB; i++) val[i] = Z_lval(N, BAD[i]);
    2047         854 :       setlg(an0, lg(vecKi)); av2 = avma;
    2048     2305926 :       for (iM = 1; iM < lM; iM++, set_avma(av2))
    2049             :       {
    2050     2305170 :         GEN M = gel(vM, iM), an, eno, den;
    2051     3214589 :         for (i = 1; i < lB; i++)
    2052             :         {
    2053     2986858 :           GEN F = gmael(vF, i, M[i]+1);
    2054     2986858 :           long D, dF = degpol(F), a = val[i];
    2055     2986858 :           if (a < d - dF) break;
    2056     2797417 :           if (a)
    2057             :           {
    2058     2606261 :             if (d == dF) break;
    2059     1163771 :             D = mael(vD, i, M[i]+1); /* sum_j deg Fj nj */
    2060     1163771 :             if (D == d && d - dF != a) break;
    2061             :             /* a = a(chi) + d - deg F - 1 */
    2062     1100575 :             if (D == d-1 && a != d - dF + achi[i] - 1 && !gequal1(t)) break;
    2063             :           }
    2064      909419 :           gmael(PPOL, i, 2) = F;
    2065             :         }
    2066     2305170 :         if (i < lB) continue;
    2067      227731 :         an = hgmmodif(an0, PPOL);
    2068      227731 :         den = gmul(t0k, RgV_dotproduct(vecK, an));
    2069      227731 :         eno = gdiv(RgV_dotproduct(vecKi, an), den);
    2070      227731 :         if (gexpo(gsubgs(gabs(eno, prec), 1)) < -bitprec / 2)
    2071             :         {
    2072          98 :           if (pole)
    2073          77 :             for (i = 1; i < lB; i++)
    2074             :             {
    2075          42 :               GEN t = deg1pol_shallow(negi(powuu(BAD[i], k2)), gen_1, 0);
    2076          42 :               gmael(PPOL, i, 2) = gmul(gmael(PPOL, i, 2), t);
    2077             :             }
    2078         231 :           for (i = 1; i < lB; i++) gmael(PPOL, i, 2) = ginv(gmael(PPOL, i, 2));
    2079          98 :           eno = grndtoi(eno, NULL);
    2080          98 :           if (typ(eno) != t_INT) pari_err_BUG("lfunhgmwild");
    2081          98 :           return mkvec3(N, mkvec2(t, PPOL), eno);
    2082             :         }
    2083             :       }
    2084         756 :       set_avma(av0);
    2085             :     }
    2086             :   }
    2087             : }
    2088             : 
    2089             : static GEN
    2090          91 : BAD2small(GEN BAD)
    2091             : {
    2092          91 :   long i, l = lg(BAD);
    2093          91 :   GEN v = cgetg(l, t_VECSMALL);
    2094         210 :   for (i = 1; i < l; i++) v[i] = itou( gmael(BAD,i,1) );
    2095          91 :   return v;
    2096             : }
    2097             : 
    2098             : /* moments of motive */
    2099             : static GEN
    2100          91 : hgmmoments(GEN H, GEN t, GEN M, long nb)
    2101             : {
    2102          91 :   pari_sp av = avma, av2;
    2103          91 :   GEN P, v, S, L = hgmlfuninfty(H, t), BAD = BAD2small(gel(L, 3));
    2104          91 :   GEN k2 = gmul2n(gsubgs(gel(L,2), 1), -1);
    2105          91 :   long ct = 0, i, j, lP, lm, tm = typ(M);
    2106             : 
    2107          91 :   if (!nb) nb = 1000;
    2108          91 :   P = primes_zv(nb); lP = lg(P);
    2109          91 :   v = hgmcoefs(H, t, P[lP - 1]);
    2110          91 :   if (tm != t_VEC) M = gtovec(M);
    2111          91 :   av2 = avma; lm = lg(M);
    2112          91 :   S = const_vec(lm - 1, real_0(DEFAULTPREC));
    2113        3731 :   for (i = 1; i < lP; i++)
    2114             :   {
    2115        3640 :     long p = P[i];
    2116        3640 :     if (!Q_lval(t, p) && !zv_search(BAD, p))
    2117             :     {
    2118        3500 :       GEN T = gdiv(gel(v, p), gpow(utoipos(p), k2, DEFAULTPREC));
    2119        3500 :       ct++;
    2120        7000 :       for (j = 1; j < lm; j++)
    2121        3500 :         gel(S,j) = gadd(gel(S,j), gpow(T, gel(M,j), DEFAULTPREC));
    2122             :     }
    2123        3640 :     if ((i & 0xf) == 0) S = gerepilecopy(av2, S);
    2124             :   }
    2125          91 :   if (tm != t_VEC && tm != t_COL && tm != t_VECSMALL) S = gel(S, 1);
    2126          91 :   return gerepileupto(av, gdivgu(S, ct));
    2127             : }
    2128             : 
    2129             : /* Heuristic guess: is there a pole ? */
    2130             : static long
    2131          91 : lfunhgmispole(GEN H, GEN t, long nb)
    2132             : {
    2133          91 :   pari_sp av = avma;
    2134             :   GEN M;
    2135          91 :   if (!nb) nb = 40;
    2136          91 :   M = hgmmoments(H, t, gen_1, nb);
    2137          91 :   return gc_bool(av, gexpo(M) > -2);
    2138             : }
    2139             : 
    2140             : static GEN
    2141         112 : tag(GEN x, long t) { return mkvec2(mkvecsmall(t), x); }
    2142             : 
    2143             : static GEN
    2144         112 : lfunhgm_i(GEN hgm, GEN t, GEN hint, long bitprec)
    2145             : {
    2146         112 :   GEN L, vr, v = hgmlfuninfty(hgm, t), vga = zv_to_ZV(gel(v,1)), k = gel(v,2);
    2147         112 :   GEN BAD = gel(v,3), COND = gel(v, 4);
    2148         112 :   long pole = mpodd(k) && lfunhgmispole(hgm, t, 0), lB = lg(BAD);
    2149             : 
    2150         112 :   L = mkvecn(pole? 7: 6, tag(mkvec2(hgm,t), t_LFUN_HGM), gen_0, vga, k, COND,
    2151             :              gen_0, mkvec(mkvec2(shifti(addiu(k,1),-1), gen_0)));
    2152         112 :   if (pole && ldata_get_degree(L) == 1)
    2153             :   { /* motive = zeta */
    2154             :     long i;
    2155           0 :     gel(L,5) = gel(L,6) = gel(L,7) = gen_1;
    2156           0 :     if (lB != 1)
    2157             :     {
    2158           0 :       GEN R = mkrfrac(gen_1, deg1pol_shallow(gen_m1,gen_1,0));
    2159           0 :       gmael3(L, 1, 2, 2) = mkvec2(t, BAD);
    2160           0 :       for (i = 1; i < lB; i++) gmael(BAD, i, 2) = R;
    2161             :     }
    2162           0 :     return L;
    2163             :   }
    2164         112 :   if (lB == 1)
    2165             :   {
    2166          14 :     vr = lfunrootres(L, bitprec);
    2167          14 :     gel(L, 6) = gel(vr, 3);
    2168          14 :     if (pole) gel(L, 7) = gel(vr, 1);
    2169          14 :     return L;
    2170             :   }
    2171          98 :   v = lfunhgmwild(L, hgm, t, BAD, pole, hint, bitprec);
    2172          98 :   gel(L, 5) = gel(v, 1); /* N */
    2173          98 :   gmael3(L, 1, 2, 2) = gel(v, 2); /* [t, PPOL] */
    2174          98 :   gel(L, 6) = gel(v, 3); /* w */
    2175          98 :   if (pole)
    2176             :   {
    2177          35 :     vr = lfunrootres(L, bitprec);
    2178          35 :     gel(L, 7) = gel(vr, 1);
    2179             :   }
    2180          98 :   return L;
    2181             : }
    2182             : GEN
    2183         119 : lfunhgm(GEN hgm, GEN t, GEN hint, long bit)
    2184             : {
    2185         119 :   pari_sp av = avma;
    2186         119 :   if (!checkhgm(hgm)) pari_err_TYPE("lfunhgm", hgm);
    2187         112 :   return gerepilecopy(av, lfunhgm_i(hgm, t, hint, bit));
    2188             : }
    2189             : 
    2190             : GEN
    2191        3997 : dirhgm_worker(GEN P, ulong X, GEN hgm, GEN t)
    2192             : {
    2193        3997 :   pari_sp av = avma;
    2194        3997 :   long i, l = lg(P);
    2195        3997 :   GEN W = cgetg(l, t_VEC);
    2196       40818 :   for(i = 1; i < l; i++)
    2197             :   {
    2198       36821 :     ulong p = uel(P,i);
    2199       36821 :     long e, d = ulogint(X, p) + 1; /* minimal d such that p^d > X */
    2200       36819 :     gel(W,i) = RgXn_inv(hgmeulerfactorlimit(hgm, t, p, d-1, 0, &e), d);
    2201             :   }
    2202        3997 :   return gerepilecopy(av, mkvec2(P,W));
    2203             : }
    2204             : 
    2205             : GEN
    2206         399 : hgmcoefs(GEN hgm, GEN t, long n)
    2207             : {
    2208         399 :   GEN worker, bad = NULL;
    2209         399 :   if (!checkhgm(hgm)) pari_err_TYPE("hgmcoefs", hgm);
    2210         399 :   if (typ(t) == t_VEC && lg(t) == 3) { bad = gel(t,2); t = gel(t,1); }
    2211         399 :   if (!is_rational_t(typ(t))) pari_err_TYPE("hgmcoefs",t);
    2212         392 :   worker = snm_closure(is_entry("_dirhgm_worker"), mkvec2(hgm, t));
    2213         392 :   return pardireuler(worker, gen_2, stoi(n), NULL, bad);
    2214             : }

Generated by: LCOV version 1.13